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ABSTRACT 


A numerical method is presented for analyzing the transonic 
potential flow past a lifting, swept wing. A finite-difference 
approximation to the full potential equation is solved in a coordi- 
nate system which is nearly conformally mapped from the physical 
space in planes parallel to the symmetry plane, and reduces the 
wing surface to a portion of one boundary of the computational grid. 
A coordinate invariant, rotated difference scheme is used, and the 
difference equations are solved by relaxation. The method is 
capable of treating wings of arbitrary planform and dihedral , 
although approximations in treating the tips and vortex sheet make 
its accuracy suspect for wings of small aspect ratio. Comparisons 
of calculated results with experimental data are shown for examples 
of both conventional and supercritical transport wings . Agreement 
is quite good for both types, but it was found necessary to account 
for the displacement effect of the boundary layer for the super- 
critical wing, presumably because of its greater sensitivity to 
changes in effective geometry. 
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INTRODUCTION 


The development of profile shapes capable of efficient 
operation in the transonic regime has spurred interest in flight 
vehicles designed specifically to operate at near sonic speeds. 

The ability to predict accurately the aerodynamic characteristics 
of the complete three-dimensional wing should have a substantial 
impact on the design of such vehicles by allowing detailed trade- 
off studies to be performed without recourse to wind tunnel test- 
ing of every design variation. 

Recent advances in the theoretical prediction of inviscid 
transonic flow fields are based largely on type-dependent, finite- 
difference solutions of the steady potential equation. These 
methods were first applied to the transonic small disturbance 
equation by Murman and Cole [1] , and the full potential equation 
by Jameson [2] and Garabedian and Korn [3] for the prediction of 
airfoil flow fields. The three-dimensional small disturbance 
equation has also been solved for swept wings by Ballhaus and 
Bailey [4] and for wing -cylinder combinations by Bailey and 
Ballhaus [5] . Finally, the full potential equation has been 
solved by Jameson for the transonic flow over an oblique yawed 
wing [6]. Although an oblique wing should be aerodynamically more 
efficient than a conventional swept wing [7] , it presents 
problems of stability and control and aeroelastic divergence. 

We consider here the prediction of the flow over a swept wing. 

In Jameson's treatment of the flow over oblique wings, 
the coordinate system is aligned in planes normal to the wing 
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leading edge. Thus, for nonzero angles of yaw the free stream 
velocity vector is not contained in these planes , and the treat- 
ment of a symmetry plane in the flow past a swept wing would be 
difficult in this coordinate system. In the analysis presented 
here, the flow is analyzed in coordinate planes parallel to the 
free stream velocity vector, and the symmetry condition is 
applied on a single coordinate surface. To allow the use of a 
fine mesh to resolve the details of the flow in the sensitive 
region near the leading edge, the spanwise coordinate lines are 
aligned with the leading edge. Thus for wings of appreciable 
sweep, the resulting coordinate system is highly nonorthogonal . 

The type of geometry we shall treat is illustrated in 
Figure 1. It consists of a wing of arbitrary planform and 
dihedral extending from a symmetry plane (or wall) . We shall 
solve a finite difference approximation to the full potential 
equation for the transonic flow past such a configuration using 
a generalized relaxation method. The finite difference approxi- 
mation is the rotated difference scheme introduced by Jameson [6] , 
and is not in conservation form. This can introduce substantial 
errors in the treatment of flows containing strong shock waves. 

To assure the correct shock jump relations one ought either to 
introduce a shock fitting scheme or else to use a difference 
scheme in conservation form. A conservative formulation of the 
small disturbance equation has been given by Murman [8] , and the 
exact potential flow equation has been solved in conservation form 
by Jameson [9] for flows past airfoils. Comparisons with 
experimental data show no clear cut advantage to using the 


conservation form without a detailed modeling of the shock wave 
boundary layer interaction [10] . This is apparently because the 
error in the shock jump relations which results from the use of 
the nonconservative schemes is in the same sense as the effect 
of the boundary layer interaction. A three dimensional scheme 
in conservation form will be discussed in a later report. 
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ANALYSIS 


Geometry 

Accurate representation of the finite difference boundary- 
conditions is much simplified if the boundairy surfaces lie in 
coordinate planes. This is achieved in the present analysis by 
a sequence of transformeitions based upon a nearly conformal mapping 
of the physical space in planes containing the wing sections , 
taken in the streamwise direction. We begin by considering the 
physical space to be described in a Cartesian coordinate system 
for which x, y, and z represent the streamwise, vertical, and 
spanwise directions, as shown in Figure 1. We then introduce an 
arbitrary singular line, just inside the leading edge of the 
profile at each spanwise station. This singular line will be the 
locus of branch points in subsequent transformations in each of 
the spanwise planes to unwrap the wing surface to a shallow bump; 
its location will be chosen to make the bump as smooth as possible. 
Representing the singular line as 

X = Xg(z) 

y = Y^iz) 

we define 

X = X - X (z) , 

o 

Y = y - yg(z) , (1) 

z = z . 

This transformation shears out the wing sweep and dihedral, and 
puts the singular line at the origin of each x,y plane. In each 
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of these planes we introduce the conformal mapping 

(Xj_ + iY^)^ = 2(x + iy) , (2) 

which maps the entire wing surface to a shallow bump near the 
plane = 0 . If we define the height of this bximp as 

Y^ = S(X,i) , 

then the final shearing transformation 

X = X^ , 

Y = Y^ - S(X,i) , (3) 

Z = z , 

reduces the wing surface to a portion of the plane Y = 0 . 

To render the computational domain finite, stretching 
transformations are introduced. For example, 

Y = — , 0. < a < 1 , (4) 

(1 - Y"^)^ 

is used to map the planes Y=+<»toY=+l. Similar transfor- 
mations are used outboard of the wing tip in the Z direction, 
and downstream of the trailing edge in the X direction. A 
sketch of the resulting rectangular computational domain is 
shown in Figure 2. 

To avoid discontinuities at the vring trailing edge, the 
branch cut in each spanwise plane is continued smoothly down- 
stream. In the physical plane, the continuation is represented 
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where t is the mean of the upper and lower surface slopes at the 
trailing edge, ^te' ^te trailing edge coordinates, and 

X is a suitably chosen scaling constant (usually taken as the 
ordinate of the local quarter-chord point) . In the solution, 
this cut is taken as the location of the vortex sheet, across 
which special difference formulas must be applied. Thus we make 
the approximation that the vortex sheet lies in a fixed surface 
near the plane of the wing which leaves the trailing edge 
smoothly according to the above formula. 


Equation of Motion 

In the absence of strong shock waves, the steady, inviscid 
motion of a compressible fluid is well approximated by the well 
known equation for the velocity potential 

(a^-u^) $ + (a^-v^) $ + (a^-w^)$ - 2uv$ - 2uw$ - 2vw$ = 0, (6) 

^ XX yy zz xy xz yz ^ ' 

where u, v, and w are the velocity components (i.e., the 
derivatives of f) in the x, y, and z directions, and a is the 
speed of sound. For the steady, potential flow of a perfect gas 
with specfic heat ratio y, '' 

2 2 Y-1 / 2^ 2. 

a = - 2“ (u + V + V7 ) , (7) 

where a^ is the stagnation speed of sound. If the flow is uniform 
at infinity, parallel to the x-y plane, and inclined at an angle a 
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to the x-axis , the far field singularity can be removed by defin- 
ing the reduced potential G as 


G = $ 

= 

The transformations of equations (1) , {2) , and (3) applied to 

equation (6) then result in an equation of the form 


- X cos a - y sin a 




r cos 




I sin a. 


( 8 ) 


If we introduce the notation 






II 

I 

X- x' - 










X 

y 



(10) 





n = 

X, x' - 
1- s 

^1. 

^s ' 








y 

X 




u = 

1 

h 

$ 

1 J 

h 1 

^X^cos 

a + Yj^s.in 

a + 

°x" ^x°y} 

t 


V = 

1 

h 

$ = 
’'i 

1 J 

h 1 

l^-Y^COS 

a+ X^sin 

a + 

s} 

r 

(11) 


I I 

w = $„ = h?U + hnV + X cos a+ y^sin a+ G„- S G , 

2 S S Zt Zi X 


and 


■/here 


U = U + hCw , 
V = V + hnw , 



d(x + iy) 
d(X^ + iY^) 


2 2 
Xi t Y^ 


( 12 ) 


(13) 
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then the coefficients in equation (9) can be written as 


2 r 2 2 

A = a + h^'? 


}- 


u 2 




B = -{a^(l+h^5^)- + (a^d+h^n^)- 

X I J 


2 2 
+ h'^Ca'^ 


-w^)S2 - |2h^a^Cn - 2Uv|s 


X 


2 ^ 2 


2 2 


+ -{2h"Ca"-2hwu}s^Sj^ - |2h‘'na‘'- 2hwV^S 


1 ^.. j^X'^Z 

' 

= - 2{ 


Z ' 


C = h 


D 


a^(l+h^?^)-U^|s^+ I 


2h^5na^- 2Uv|- |2h^5a^“ 2hwu|s2 , 


E = 2h^Ca^ - 2hwU , 


{a^(l+h^5 ^)- U^l^xx" h^(a^-w^)S 


|2h^5a^-2hwu|s_.- >-G 


XZ Y 


F = -2h^ (a^-w^)S 2 - |2h^Ca^- 2hwu|s^+ 2h^na^- 2hwV , 

R = {-{a 

+ h 3 (a 2 _„ 2 ,f|(^. 2 .y. 2 , 2 x>;x - x^X - y”Xj_}u 

'> XX xy X 

+ 2 x;y^x + x;x - y^x^ }v} 

xy XX y X"' ^ 

+ 2h^((X x;-X y;)X + (X xV X^_y;)Xj^_.}(U^+v2) 

*' X y XX y X xy-* 

■*■ E {x^U+y^v|(U^+V^) + cos a|h^(e^-n^)a^-U^-V^+h^(a^-w^)x^| 
in 


(14) 


+ sin a-{2h^Cna^- 2UV + h^(a^-w^)yg| . 


Note that for the transformation defined by equation (2) , 


X^_ = X^/h' 
Y^_ = Y^/h‘ 

y 


(15) 
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and 


^1~ 

XX 


12 2 

I (h^ - 4Y^) , 


(16) 


X 


^xy 


I (h^ - 4X^) . 


The symmetry condition that w = 0 on the plane z 
requires 


= 0 




- {"z - 'll 


Gy = 0 , 


(17) 


and the boundary condition that the flow be tangent to the wing 
surface requires 


h + ®X> + {®Z ^ SSj, 

n 



+ 


+ 


{ ^2 ^x'*’ ^^X‘*'^}}^X { 

■I -X- _ COS a- X^_sin ot j-S - X,_cos 

^ ■^x y y 


=z- i}s 

a+X^_sin a - 0 , 
■^x 


(18) 


on Y = 0. 

Downstream of a finite lifting wing there will be a vortex 
sheet. Across the sheet the pressure is continuous, but there 
may be discontinuities in the tangential velocity components. 
Convection and roll-up of the vortex sheet are ignored. In 
reality, the component of velocity normal to the sheet must be 
zero, but in our approximation it is simply required to be contin- 
uous . Thus , the equation 

♦yy = ° 

is used at points lying on the vortex sheet. Also the disconti- 
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nuity in potential is assumed to be constant along streamwise 
coordinate lines downstream of the trailing edge. The value of 
this discontinuity is determined by the Kutta condition^ and its 
spanwise variation determines the strength of the vortex sheet. 

Finite Difference Approximation 

The success of the type dependent difference scheme applied 
to the transonic small disturbance equation by Murman and Cole [1] 
can be attributed to the fact that it effectively adds a direc- 
tional bias to the equation at points where the local flow is 
supersonic. In constructing an analogous scheme for the full 
potential equation in general curvilinear coordinates (which may 
not be aligned, even approximately, with the local flow direction), 
care must be taken to ensure that this bias is added in the upwind 
direction, i.e., in the direction parallel to the velocity vector. 

A method with this property has been proposed by Jameson [6] . 
To illustrate it, we return to the potential equation in the 
physical coordinates. The equation is rearranged as if it were 
expressed in a Cartesian coordinate system aligned with the local 
flow direction, s, at the point under consideration. Then 
equation (6) assiames the canonic?.! form 

(a^-q^)$gg + a^(V^$-$gg) = 0 (19) 

where q is the magnitude of the velocity. 
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The relaxation scheme is designed to simulate an artificial 
time dependent process which converges to the desired solution 
of the steady state equation. In the finite difference approxima- 
tion to the potential equation, central differences are used to 
calculate all first derivatives, from which the velocities can be 
determined using equations (11) . At grid points where the flow 
is subsonic, central differences are also used to approximate the 
second-order derivatives in equation (9) . A typical central 

difference formula for is 

jCX 

(n+1) 2 (n+l) (n) (n) 

G = ^i-l,j,k ~ ^i,j,k ~ *^i,j,k ^i-H,j,k 

■ ax2 


( 20 ) 


Where the superscripts denote the iteration level and w is 
the relaxation factor [6] . If we regard each iteration as repre- 
senting an advance At in an artificial time coordinate, this 
formula can be interpreted as an approximation to 


XX 


if ^^Xt Ax ^03 


Similarly, the formula 

^(n) p(n) p(n+l) „ (n+1) 

r - i-tl,j+l.k ~~ ^i+l,i-l,k ^i-l,j+l,k ^ ^i-l.j-3,k 

^XY 4AXAY 


(21) 


can be interpreted as an approximation to 

1 At ^ 


XY 


2 Ax ^Yt * 


The relaxation process can thus be regarded as an approximation to 
the time dependent equation 


11 


(M -1)G - G - G + 2a,G 2a^G 2a.G .+ 6G^ 

ss Iran nn 1 st 2 mt 3 nt t 


where M = q/a is the local Mach number, m and n are suitably ; 
scaled coordinates in the plane normal to the velocity vector,; 
and Q contains all the terms in the equation other than the j 
principal part. The coefficients a^,a 2 / 0 t 3 / and 5, depend on ; 
the mix of old and updated values in the difference equations 
as well as any explicit time-like or mixed terms that have be(Sn 
added for stability. 

Introducing the new time coordinate 

“l 

T = t - —2 — s + tt2m + ' 

M —1 

transforms equation (22) to 


► r *^1 2 2l 

'-DG - G - G - 

ss mm nn |^j^2 ^ 2 3J 


^TT ^ ^ * 


In order to ensure the convergence of the scheme, we require 
that equation (23) should be a damped three-dimensional wave 
equation. This will be the case if 

> (M^-1) 

At points where the velocity is supersonic, upwind differ- 
ences are used to represent contributions to G in the first 
term of equation (19) . This is done using formulas of the type 


2Gf^t^^ - Gf^^. , 

i,3,k i,n,k 


(n+1) . p (n) 

^^i-l,j,k ^ ^i-2,j,k 


p(n+l) p(n+l) p(n+l) . p (n+1) 

^i,i,k ~ ^i-l,j,k ~ ^i,i-l,k ^ ^i-l,j-l,k 

AXAY 
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_fc. 


JU. 




These formulas also have the property of guaranteeing diagonal 
dominance for the updated values on each line. The formula for 
G can be interpreted as representing 


^XX ^ AX ^Xt 


Together with analogous formulas for and G 

X X h L 

duces a term, equal to 




this intro- 


into equation (22) . To ensure that equation (24) is satisfied 

2 

at points near the sonic line where (M -1) is small, the coeffi 

cient of G , can be further augmented by adding a term of the 
s ^ 

form 

B n Ht ^ ^ ' 

where 3 > 0 is appropriately chosen. The required mixed 

derivatives can be constructed in the form 

(n+1) _ (n) _ p(n+l) ^,(n) 

At p ^ i,j,k i,j,k i-l,j,k i-l,j,k 
AX >^Xt 


The supersonic difference scheme is completed by using central 
difference formulas similar to equations (20) and (21) to 
evaluate contributions to the second term of equation (19) , 
but with 0 ) set to unity, as suggested by a local von Neumann 
test [6] . 


I 

* 


(26) 


(27) 
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Boundary Conditions 


The boundary condition at infinity is particularly simple 
because the square root transformation reduces the entire vortex 
wake to the X-Z plane at downstream infinity. Therefore, since 
the uniform stream singularity has been removed by the introduc- 
tion of the reduced potential, the Dirichlet condition 

G = 0 

is appropriate. 

On the X-Y and X-Z planes, finite difference approximations 
to the Neumann boundary conditions specified by equations (17) and 
(18) must be applied to those portions representing solid bounda- 
ries (i.e. , the symmetry plane and the wing surface) . At the wing 
surface, central difference approximations are used in equation 
(18) to define values of the reduced potential at image points 
located one mesh spacing below the X-Z plane. A similar method 
is used on the symmetry plane, but due to the high degree of 
nonorthogonality of the coordinate system when the wing is highly 
swept, simple central differences become unstable. Thus, to set 
the potential values at the image points for the symmetry plane, 
the X-dif f erences required in equation (17) are evaluated by 
averaging one-sided differences on either side of the symmetry 
plane, taken in the upwind direction in the image plane, and in 
the downwind direction in the first plane in the flow region. 

The symmetry condition thus remains formally second order accurate, 
and the incorporation of the image point whose value is being set 
into the X-dif ference adds to the stability of the scheme. This 
method of handling the symmetry condition has proved stable for 
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sweep angles in excess of 35 degrees. 

At points on the X-Z plane which do not lie on the wing 
surface, the values of the reduced potential at the image points 
are taken to be those of the associated point on the other side 
of the branch cut, allowing for a discontinuity across the vortex 
sheet. The value of this discontinuity is taken to be independent 
of X at each spanwise station, and its value is determined by the 
Kutta condition that the flow leave the trailing edge smoothly. 

One final note concerns points which lie on the contin- 
uation of the singular line outboard of the wing tip. At 
these points the mapping is singular, and a special 
limiting form of the difference equations must be used. At 
points where the solution is regular, the nonlinear terms of 
the potential equation are of 0(l/h), while the Laplacian 
transforms to 


, 2 ^ X,X. -ZZ 

h 11 11 

Thus, in the limit as h tends to zero. 


$ + * 1 = 0 


(28) 


is a suitable limiting form. 
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RESULTS 


Computational procedure 

The potential formulation is particularly attractive for 
three-dimensional calculations because it requires the storage of 
only one quantity at each grid point, and the number of grid 
points required to accurately describe these flow fields is large 
Even so, it is impractical to store the entire solution array 
in the high speed core of many current computing machines . 
Fortunately, since the analysis presented here depends on a 
relaxation solution of the difference equations, it is not neces- 
sary to have the entire solution immediately available at all 
times. It is, therefore, stored on a disk file, and read into 
core one X-Y plane at a time. At any time during the solution 
procedure, the values of the potential on four such planes are 
in the core. Old values are buffered in and new values buffered 
out of core while other calculations are being performed as much 
as possible, to keep the process efficient. 

In each X-Y plane, the equations are solved by successive 
line overrelaxation. The plane is divided into three regions, 
as shown in Figure 3. In the central region the equations are 
relaxed along horizontal lines , sweeping from infinity to the 
wing surface. In the outer regions the equations are relaxed 
along vertical lines, sweeping away from the central region to 
infinity. Such a sweep pattern ensures that the sweep direction 
will not be opposed to the flow direction in any supersonic zones 


16 



which would result in instability. In many cases, the central 
region can be taken to cover the entire plane; that is, only 
horizontal line relaxation is used. 

To speed convergence, an initial calculation is usually 
performed on a coarse grid, typically containing 48x6x8 grid cells 
in the X, Y, and Z directions respectively. This solution is 
then interpolated onto a finer grid containing twice as many mesh 
cells in each direction, and is used as a starting guess for an 
intermediate solution. The process is repeated once again to 
give the final solution on a grid containing 192x24x32 mesh cells. 
A typical run consists of 100 relaxation sweeps on each grid, 
requiring a total of approximately 85 minutes of CPU time on a 
CDC 6600. The same program has been run on the CDC 7600, for 
which a similar calculation requires about 15 minutes. 


Examples 

In this section we present the results of calculations 
using the swept wing program, and compare the predicted surface 
pressure distributions with those measured in experiments. The 
comparisons are made for two different wings, each typical of 
a class of swept wings of the subsonic transport type. 

The first wing geometry is representative of the tip 
panel of a relatively simple wing of conventional high speed sec- 
tion shape. It has a uniform section of 9.8 percent thickness ratio. 
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and the planform has a leading edge sweep angle of 30°, a taper 
ratio of 0.7, and an aspect ratio of 3.8. A program generated 
projection drawing of the wing is shown in Figure 4. The wing 
was tested by Monnerie and Charpin [11] of the ONERA, and carries 
their designation of wing M-6. 

The first results presented are at a free stream Mach 
nixmber of 0.9226 and zero angle of attack, resulting in zero lift 
for this symmetrical wing. Figure 5 compares the calculated and 
measured streamwise surface pressure distributions at the 20, 45, 
65, and 95 percent semispan locations [11,12]. Agreement is 
quite good, including the predxcted shock location. 

Figure 6 shows similar results for the same wing at a 
Mach number of 0.919 and an angle of attack of 3.07 degrees. 

Again, agreement between the computed and experimental results 
is quite good, with the exception of the shock location on the 
lower surface, which is somewhat further aft than predicted by 
the calculation. 

Figure 7 shows a program generated, three-dimensional, 
projection view of the wing surface pressure distribution at a 
Mach nximber of 0.840 and an angle of attack of 3.06 degrees. 

This is a particularly interesting case because of the merging 
of two shocks into one on the wing upper surface as one proceeds 
outboard. This pattern is graphically illustrated in the projec- 
tion view. Figure 8 shows comparisons of the calculated results 
with experiro,ental data, again at the 20, 45, 65, and 95 percent 
semispan stations. Agreement is quite good, including the 
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prediction of the double-shock pattern at the inboard stations . 

Figure 9 shows the projection view of the wing surface 
pressure distribution at a Mach number of 0.837 and an angle 
of attack of 6.06 degrees. Again, the calculation predicts the 
merging of a double shock pattern inboard to a single shock 
further outboard. Comparisons with data, shown in Figure 10 
show that agreement is still quite good. 

The second geometry is representative of wings being 
considered for the next generation of subsonic transport aircraft. 
The wing is twisted, both aerodynamically and geometrically, is 
highly tapered, and has a discontinuity in trailing edge sweep 
angle at the 35 percent semispan location. The planform has a 
leading edge sweep angle of 35 degrees and an aspect ratio of 7. 

It has 5 degrees of dihedral. It is defined by four distinct 
streamwise sections (at the 12, 35, 70, and 100 percent semispan 
stations) , with linearly interpolated coordinates between. The 
streamwise thickness ratio varies from 16.3 percent at the root 
to 11.9 percent at the tip. For the wind tunnel tests the wing 
was mounted on a quasicylindrical fuselage which extended to 
the 12 percent semispan. For the computations, the symmetry 
plane was assumed to be at the same spanwise station as the 
wing-fuselage intersection in the tests. A projection drawing 
of the wing (extended to the fuselage centerline) is shown in 
Figure 11. For these calculations, the wing geometry was modified 
to account for boundary layer effects by adding the displacement 
thickness obtained from two-dimensional boundary layer calculations 
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multiplied by an empirically determined spanwise weighting 
factor. The wing was one of several tested in a cooperative 
program by the Douglas Aircraft Company and the NASA Ames 
Research Center in the Ames 11-foot tunnel at a Reynolds number 

g 

of approximately 5xlo , based on the mean aerodynamic chord. 

A program generated three-dimensional projection drawing 
of the upper and lower surface pressure distributions for this 
wing is shown in Figure 12. (This particular case was . run with 
no correction for boundary layer displacement effect, and with 
the wing extended to the fuselage centerline.) 

Comparisons with experimental data are shown in Figures 
13 and 14. The first case. Figure 13, shows streamwise surface 
pressure distributions at a number of spanwise stations for a 
Mach nximber of 0.75 and an angle of attack of 2.2 degrees. 
Agreement with experiment is seen to be excellent, including 
the location and strength of the rather strong shock near the 
leading edge on the wing upper surface. 

Figure 14 shows similar comparisons at a Mach number of 
0.84 and an angle of attack of 1.85 degrees. Again, agreement 
is quite good, although the resolution of the first (rather weak) 
shock of the inboard double shock pattern seems lost between the 
35.5 and 50 percent semispan locations. 

The results displayed in Figures 13 and 14 were kindly 
supplied by R. M. Hicks and P. A. Henne. Further details of 
the wing geometry, calculations, and test conditions are 
contained in [13] . 
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CONCLUSIONS 


A numerical method has been presented for determining the 
inviscid transonic flow past a swept wing. The method is based 
on a type-dependent, finite difference approximation to the full 
potential equation, solved in a computational domain designed 
for accurate application of the wing surface and symmetry plane 
boundary conditions. Calculated surface pressure distributions 
agree well with experimental data for wings of conventional and 
supercritical section shape (when the geometry in the latter 
cases is corrected for the displacement effect of the boundary 
layer) . 

Mapping techniques similar to those used here could be 
used to treat more realistic geometries, e.g., a wing mounted 
on a fuselage [14] . The recasting of the finite difference 
approximation into conservation form would also be an important 
theoretical contribution. 

Finally, as was mentioned in the preceding section, these 
calculations require a substantial amount of computer time. 

Thus , methods of accelerating the convergence of the iterative 
scheme are particularly important in three-dimensional problems. 

A number of techniques to achieve this have met with success in 
two-dimensional calculations, including a hybrid Poisson-solver/ 
relaxation technique [15,16] , a multi-grid method [17] , and an 
alternating-drection method [18] . The extension of these methods 
to three-dimensional calculations should result in great savings. 
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(a) Plan View 


(b) Front View 


Figure 1. Geometry of Swept Wing, 
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Figure 5. Comparison Of Calculated And Experimental Wing Pressure Distributions For Onera Wing M"6 
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Figure 6. Comporison Of Calculated And Experimental Wing Pressure Distributions For Onera Wing M"6. 





FIGURE 7 


UPPER SURFfiCE PRESSURE LOWER SURFfiCE PRESSURE 

ONERfl WING M6 L* E. SWEEP 30. DEG. fiSPECT RfiTIO 3.8 
MflCH .840 YRW 0.000 RLPHR 3.060 

L/0 13.89 CL .2860 CD .0206 
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FIGURE 9 


UPPER SURFACE PRESSURE LOWER SURFACE PRESSURE 

ONERA WING M6 L. E. SWEEP 30. DEG. ASPECT RATIO 3.8 
MACH .837 YAW 0.000 ALPHA 6.060 

L/D 9.61 CL .5587 CD .0581 
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FIGURE 12, THREE-DIMENSIONAL SURFACE PRESSURE DISTRIBUTION. 

UPPER SURFACE PRESSURE LOWER SURFACE PRESSURE 

DOUGLAS WING W2 (EXTENDED TO CENTER LINE) 

MACH .819 YAW 0.000 ALPHA Q.OQO 

L/D SO. 09 CL .5155 CD .0272 








Figure 14. Comparison Of Calculated And Experimental Wing Pressure Distributions For DAC Case 5 






Appendix A. Description of the program 


All the numerical results in this report were generated 
by the computer program FLO 22 listed in Appendix B. This 
program includes options to treat both a swept wing on a 
wall (Figure Al) , and an isolated yawed wing (Figure A2) . 

For swept wing calculations the sheared parabolic coordi- 
nates are introduced in planes parallel to the free stream. 

In the treatment of a yawed wing the whole coordinate system 
is rotated through a specified yaw angle, so that the X-Y 
planes are normal to the leading edge of the wing at its 
center line. In either case the wing section can be varied 
in an arbitrary manner, and the only restriction on the plan- 
form is that the leading edge may be any smooth curve, but 
it should not have kinks, since these would cause the second 
derivatives of the singular line of the coordinate system 
to become unbounded. Kinks are permitted in the trailing edge, 
on the other hand. The trailing edge defined by the input 
is actually replaced by a piecewise straight line connecting 
the nearest mesh points in the computational lattice . 
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The geometry is defined by giving the wing sections at 
successive span stations from the wing root to the tip, or in 
the case of a yawed wing, from the leading to the trailing tip. 

Up to 11 span stations may be used for this purpose, and the 
planform and dihedral are determined by specifying the chord 
and the x and y coordinates of the leading edge at these span 
stations. The wing section at each station is then determined 
by scaling and rotating a prescribed profile, given by a table 
of X and y coordinates. If the wing sections are similar, only 
the profile for the first station need be read in. The coordi- 
nates for the other stations are obtained by scaling the original 
profile to the proper chord, and rotating it to obtain the 
appropriate twist. If, on the other hand, the sections are not 
similar, the program permits the coordinates of new profiles to 
be read in at each span station. The wing section between 
stations is generated by interpolation. The location of the 
singular line about which the wing is unwrapped by the square 
root transformation is determined by the parameters XSING and YSING, 
which must be specified at each span station. It is important 
to choose these so that the mapped profile does not have any sharp 
bvimps . 

The main input to the program is read from Tape 5, and 
the output is written on Tape 6. Tapes 1, 2 and 3 are disk files 
used for internal storage in order to reduce the requirements for 
high speed memory. Tape 4 is a permanent storage device such as 
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a magnetic tape cn which an intermediate result can be saved. 

The computation can then be continued for more iterations , 
starting from the values saved on Tape 4. The disk instruc- 
tions in the version of the code listed in Appendix B are 
specialized to the CDC 6600 using the FTN compiler. Otherwise 
the code should be readily adaptable to other computers. 

The data deck for a run is arranged to include 
title cards listing the required data item.s . The complete set 
of title cards provides a list of all the data which must be 
supplied, and can be used as a guide in setting up a data deck. 
Each title card is followed by one or more cards supplying the 
numerical values of the parameters listed on the title card. 

All data items are read as floating point niombers in fields of 
10 columns, and values representing integer parameters are 
converted inside the program. A glossary of the input parameters 
is given in Table 1, and a typical data deck is shown in Table 2. 
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Table 1. Glossary of input parameters 

(Listed in order of their occurrence on the data title cards) 


TITLE CARD 1 
NX 

NY 

NZ 

FPLOT 

XSCAL, PSCAL 


FCONT 


The number of mesh cells in the direction of the 
chord used at the start of the calculation. 

NX = 0 causes termination of the program. 

The number of mesh cells in the direction normal 
to the chord and span. 

The number of mesh cells in the span direction. 
Controls generation of plots. 

FPLOT=0 . for a print plot but no Calcomp plot 
at each span station. 

FPL0T=1. for both a print plot and a Calcomp plot 
at each span station. 

FPL0T=2 . for a Calcomp plot but no print plot at 
each span station. 

FPL0T=3 . for a three dimensional Calcomp plot only. 

Control the scales of the Calcomp plots . 

XSCAL>0. scales each section plot to XSCAL 
XSCAL=0 . scales each section plot to 5.0 
XSCAL<0. scales the maximum chord to XSCAL, and 
each section plot proportionately to the local chord 
PSCALt^O. sets the pressure scale to PSCAL per inch 
in each section plot. 

PSCAL=0 . sets the pressure scale to 0.4 per inch 
in each section plot. Also, 

PSCAL^O. scales the three dimensional plot so 
that the span or semispan is 5. If PSCAL=0 . and 
XSCAL^^O. then the three dimensional plot is 
scaled so that the maximiam chord is 1/2 XSCAL. 

Indicator which determines the manner of starting 
the program. 

FCONT=0. indicates the calculation begins at 
iteration zero. 

FC0NT=1. indicates the computation is to be 
continued from a previous calculation. In this 
case the values of the velocity potential and the 
circulation are read from a magnetic tape where 
they were previously stored (Tape 4) . It is still 
necessary to provide the complete data deck to 
redefine the geometry. The count of the iteration 
cycles is continued from the final count of the 
previous calculation and the maximum number of 
additional iterations to be performed is defined 
by MIT. 
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TITLE CARD 2 
MIT 

COV 

PI 

P2 

P3 

BETA 

STRIP 

FHALF 


The maximum number of iteration cycles which will 
be computed. 

The desired accuracy. If the maximum correction 
is less than COV the calculation terminates or 
proceeds to a finer mesh, otherwise the number 
of cycles set by MIT are completed. 

The subsonic relaxation factor for the velocity 
potential. It is between 1. and 2. and should 
be increased towards 2. as the mesh is refined. 

The supersonic relaxation factor for the velocity 
potential. It is not greater than 1. and is 
normally set to 1. 

The relaxation factor for the circulation. 

It is usually set to 1., but can be increased. 

The damping parameter controlling the amount of 
added (see equation (2.6), page 13). 

It is normally set between 0. and 0.25. 

Determines the split between horizontal and 
vertical line relaxation and is the proportion 
of the total mesh in which horizontal line relaxa- 
tion is used. Fastest convergence is usually 
obtained by setting STRIP =1. so that horizontal 
line relaxation is used for the entire mesh. 

If convergence difficulties are encountered STRIP 
may be reduced to some fraction between 0 . and 1 . 

Determines whether the mesh will be refined. 
FHALF=0.: the computation term.inates after 
completing the prescribed number of iteration cycles 
or after convergence. 

FHALFt^O . ; the mesh spacing will be halved after MIT 
cycles have been run on the crude mesh size. An 
additional data card must be provided for the 
refined mesh giving the numerical values requested 
by Title Card 2. If 

FHALF<0 the interpolated potential will 
smoothed j FHALF [ times . 
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TITLE CARD 3 


FMACH 

YAW 

ALPHA 

CDO 

TITLE CARD 4 
ZSYM 

NC 

SWEEP 1 
SWEEP 2 

SWEEP 

DlHEDl 

DIEED2 

DX'HED 


The free stream Mach number. 

The yaw angle of the wing in degrees . 

The angle of attack in degrees . When the wing 
is yawed, ALPHA is measured in the plane normal 
to the leading edge, not in the free stream 
direction. 

The estimated parasite drag due to skin friction 
and separation. It is added to the pressure drag 
(sum of vortex drag plus wave drag) calculated 
by the program to give the total drag. 


Determines whether to treat a wing on a wall or 
an isolated wing. 

ZSiM=l.; the wing is on a wall 

ZSYM=0.; the wing is an isolated wing at a yaw 

angle given by YAW. 

The number of span stations at which the wing section 

is defined on subsequent data cards from the wing 

root to the tip if ZSYM=1., or from the leading 

to the trailing tip if ZSYM=0 . If 

NC<3 it is assumed that the wing geometry is 

the same as for the last case calculated and 

the computation for new values of FMACH, YAW, ALPHA 

and CDO begins without further data items 

being read. 

Sweep of singular line at the v/ing root if ZSYM=1., 
or at the leading tip if ZSYM=0 . 

Sweep of singular line at the tip. 

( SWEEP 1 and SWEEP 2 are used as end conditions 
for a spline fitting the x coordinates of the 
singular line.) 

Sweep of singular line in the far field. 

Dihedral of singular line at the wing root if 
ZSYM=1., or at the leading tip if ZSYM=0 . 

Dihedral of singular line at the tip. 

(DlHEDl and DIHED2 are used as end conditions for 
a spline fitting the y coordinates of the singular 
line . ) 

Dihedral of singular line in the far field. 
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TITLE CARD 5 
Z 

XLE , YLE 
CHORD 

THICK 

ALPHA 

FSEC 


TITLE CARD 6 
YSYM 


NU 

NI, 

TITLE CARD 7 
TRAIL 

SLOPT 


/ 


(The geometry at the first span station) 

Span location of the section. 

X and y coordinates of the leading edge. 

The local chord value by which the profile 
coordinates are scaled. 

Modifies the section thickness. The y coordi- 
nates are multiplied by THICK. 

The angle through which the section is rotated to 
introduce twist. In the case of a yawed wing^ this 
angle is measured in the axis system attached to 
the wing, not in the direction of the free stream. 

Indicates whether or not the geometry for a new 
profile is supplied. 

FSEC=0 . : the section is obtained by scaling 
the profile used at the previous span section 
according to the parameters CHORD, THICK, ALPHA. 

No further cards are read for this span station, 
and the next card should be the title card for the 
next span station, if any. 

FSEC=1.: the coordinates for a new profile are 
read from the data cards which follow. 

(Profile Geometry Supplied if FSEC=1.) 

Indicates the type of profile. 

YSYM=0. denotes a cambered profile. Coordinates 
are supplied for upper and lower surfaces, each 
ordered from nose to tail with the leading edge 
includtsd in both surfaces. 

YSYM=1. denotes a symmetric profile. A table 
of coordinates is read for the upper surface only. 

The mrniber of upper surface coordaintes . 

The number of lower surface coordinates. 

For YSYM=1., NL=NU even though no lower surface 
coordinates are given. 

(Additional Profile Geometry Supplied if FSEC=1.) 

The included angle at the trailing edge in degrees. 
The profile may be open, in which case it is the 
difference in angle between the upper and lower 
surfaces . 

The slope of the mean camber line at the trailing 
edge. This is used to continue the coordinate 
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-t: 


XSING, YSING 


TITLE CARD 8 
X,Y 

TITLE CARD 9 
X,Y 

TITLE CARD 10 


surface, assumed to contain the vortex sheet, 
smoothly off the trailing edge. For heavily aft 
loaded airfoils, the lift is sensitive to the 
value of this parameter, which should be adjusted 
by comparing two dimensional calculations using 
parabolic coordinates with two dimensional calcula- 
tions in the circle plane. 

The coordaintes of the singular point inside the 
nose about which the square root transformation 
is applied to generate parabolic coordinates. 

This point should be located as symmetrically as 
possible between the upper and lower surfaces at 
a distance from the nose roughly proportional to 
the leading edge radius. It can be seen whether 
the location has been correctly chosen by inspect- 
ing the coordinates of the mapped profile printed 
in the output. If the mapped profile has a bump 
at the center, the singular point should be 
moved closer to the leading edge. If the mapped 
profile is not symmetric near the center, with a 
step increase in y, say, as x increases through 0, 
the singular point should be moved closer to the 
upper surface. The coordinates of the singular 
point are chosen relative to the profile coordinates 
supplied on the cards which follow. 

(Upper Surface Coordinates) 

The coordinates of the upper surface. These are 
read on the data cards which follow, one pair of 
coordinates per card in the first two fields of 10 , 
from leading to trailing edge inclusive. 

(Lower Surface Coordinates, Read if ISYM = 0.) 

The coordinates of the lower surface, read from 
leading edge to trailing edge. The leading edge 
point is the same as the upper surface leading edge 
point. The trailing edge point may be different if 
the profile has an open tail. 

11... (Geometry at the Other Span Stations) 

These title cards are the same as Title Card 5 
(geometry for the first span station) . The number 
of such cards depends on the number of input span 
stations NC. If the profiles are similar at each 
station except for scaling, thickness to chord ratio 
and rotation to introduce twist, FSEC=0 . and no 
new profile coordinates are needed. 
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TABLE 2. DATA DECK FOR ONERA M6 WING 


Columns 

Cards 


1-10 1 11-20 


21-30 


31-40 i 41-50 I 51-60 


Title of case ONERA M6| WING (copied onto output ahd Calcom^ plots) 


Title Card 


Title Card 


Title Card 




FPLOT 

1 . 


Title 

Card 

MIT 

COV 

f==== 

PI 

P2 i 

P3 

1 

BETA i 
1 1 

STRIP 

1 

FHALF 



100. 

l.E-6 

1.6 

1. i 

1. 

‘ I 

.10 1 

1. 

1. 



100. 

l.E-6 

1.6 

1. 

1. 

.10 . 

1. 

1. 



100. 

l.E-6 

1.6 

1. 

1. 

.10 . 
« 

1. 


Title 

Card 

MACH 

YAW 

ALPHA 

CDO 


( 





.840 

0. 

3.06 

.010 





Title 

Card 

ZSYM 

NC 

SWEEP 1 

SWEEP 2 

SWEEP 

DIHEDl ! 

J 

DIHED2 

DIKED 



1. 

6. 

29.9 

29.9 

29.9 

0. i 
1 

0. 

1 

0. 

Title 

Card 

Z 

XLE 

YLE 

CHORD 

THICK 

ALPHA 

FSEC 




0. 

0. 

0. 

.6737 

l’. 

0. 

t 

1. 




YSYM 

1 . 


TRAIL SLOPT XSING YSING 



,00725 0 









XLE 

.1150 


XLE 

.2300 


XLE 

.3450 


XLE 

.4600 


XLE 

.5750 







CHORD 

.6147 


CHORD 

.5558 


CHORD 

.4968 


CHORD 

.4379 


CHORD 

.3789 


THICK 

1 . 


ALPHA 

0 . 


ALPHA 

0 . 





THICK I ALPHA 



THICK 

1 . 


i ALPHA 
! 0 . 



THICK ALPHA 
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Both graphical and printed output are provided. The 
wing sections defining the geometric configurations are 
printed for each span station/ if they are different, or for the 
first span station only if the sections are all similar. The 
program next prints the coordinates of the unfolded sections 
produced by the square root transformations at the root and 
the tip. These should be inspected to see that they are reason- 
ably smooth. The program also prints a chart of an indicator IV 
showing the configuration of the wing in the coordinate surface 
to which it has been mapped. The values of IV are as follows: 

IV = 2 indicates a point on the wing 

1 indicates a point on the trailing vortex sheet 

0 indicates a point on the singular line 

-1 indicates a point adjacent to the edge of the wing 

or vortex sheet 

-2 indicates an ordinary point not in contact with the 
wing or vortex sheet. 

The program next displays the iteration history. The 
maximum correction to the velocity potential and the maximum 
residual of the difference equations are printed at each cycle, 
together with the locations of the points where these occur 
in the computational lattice, and also the relaxation factors, 
the circulation at the wing center line, and the number of 
supersonic points. 
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After a specified maximum number of cycles has been 
completed, or a convergence criterion has been satisfied, the 
section lift, drag and moment coefficients are printed for 
each span station, and the pressure distribution is printed 
or displayed in a Calcomp plot as desired. Finally the charac- 
teristics of the complete wing are printed. These include the 

coefficients of lift and form drag computed by integrating 
the surface pressure, and the ratio of lift to form drag. 

An estimate of the friction drag coefficient may be supplied 
in the input, and this vjill be included to provide an estimate 
of the total drag coefficient of the ratio of lift to total drag. 
The pitching, rolling and yawing moments are also computed 
and printed. In the case of a yawed wing these are in an axis 

system normal to the wing leading edge at its center line. In 

the case of a wing on a wall the rolling moment is the root 
bending moment. 

Finally additional Calcomp plots are generated if they 
are desired. These show the convergence history, and also a viev7 
of the complete wing and the three dimensional pressure distri- 
bution over the upper and lower surfaces separately, with the 
wing root or the leading tip at the bottom of the picture. If 
the mesh is to be refined the program, then completes the same 
sequence of calculations and output for the new mesh. 
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APPENDIX B. LISTING OF THE PROGRAM 


PROGRAM FL 02 ? ( INPUT f OUTPUT, TAPE i# Tape 2 ^ TAP t 3 > TAP 54 ^ 

TAP Ei»» INPUT, TAP t6 -OUTPUT) 

THREE DIMENSXGNAL wIN 6 ANALYSIS IN TkANSQMC FLOw 

USING SHEARED PARAeOlIC COOkCINATES 

WITH STORAGE ON THE DiiC 

PKGGPAMHEO 8Y ANTONY vAMb SON, PARCH 197^ 

REVISIONS 6 Y D. A. CAUGHtY ANL ANTONY JA'lESON,JEC 1973 -DbC 1976 
G IS REDUCED VELOCITY POTENTIAL 


Cum. 'I ON 


2 

3 

4 

6 

7 


COMMGN/FLO/ 

dimension 


<1 

3 

4 

5 

6 
7 


ND 

Nt 

IREAO 

IWRIT 

KPLOT 

IPLQT 

ISTQP 

N 1 

N 2 

N 3 

REWIND 

REWIND 

REWIND 

REWIND 

JC 

RAO 

WRITE 

WRITE 


G(i9 3,, 2 6, ‘♦),SO(193, 35), EO( 131), ZO( 131), 

IUi9 3, 25),nt:l{35),ITt:d(35), 

AC'(193),Al(l9 3),A2(19 3),A3(19 3), 

BG( 26), 81(26), 62(26), 83(26), 

Z(35),C1( 3:,),u2(35),C3(35), 

XC(35),XZ(35),X2Z(35),YC(35),YZ(35),YZZ(35), 

NX, N Y,*^Z,KTE.l,K.Tt2, iSYM,KSYM,SCAL,SCALZ, 

YAW,CYAW, SYAi,, ALPHA,CA, SA,FMACH,n 1,N2, N3, 10 
STRIP,P1, P2,P3,BETA,FR, IR, JP,KR,OG, IG, JG,KG,NS 
XS(2A1,11), YS(2Al,il), 

ZSdl ), XLE(lx),YlE(il), SLOP 1(11), TRAIL (11), NP (II), 

El(ll),t2(li),L3(ll),F4(ll),b6(ll), 

XP(241),YP(2Al),Dl(2Ai),D2(2Al),03(241), 

X(193),Y(i93),SV(193),SM(193),CP(193), 

ChCRD(35),SCL(35),SC0(35),SC.M(35),TlTLE(20), 

FIT(3),CuVu(3),P10(3),P20(3),P30(3),B£TA0(3), 

ST^iPOl 3 ) , FHaLF( 3 ), RtS( 501),CGUNT( 601) 

2A1 

193 

5 

6 
C 
1 
2 
1 
2 
3 


1 

2 

3 

4 


„ rt-V OF 


0 


= 57.2957795130823 
{ IwRIT,6C0) 

( IWRIT, 2 ) 

FORMAT { 14H0PP0GPAM F L 02 2, 70X, 32 H AN TON Y JAMESON, COURAN T INSTITUTE/ 
50HOTHREE DIMENSIONAL WING ANALYSIS IN TRANSONIC FLOW, 

36H USING SHEARED PARABCLIC COORDINATES) 


READ 

( IPEAD, 530) 

TIILE 

WRITE 

( 1WR1T,630) 

TITLE 

READ 

( IRfc AD, 50t ) 


READ 

( I read, 510 

FNX, 

NX 

» FNX 



FNX, FNY, FNZ, FPLOT, XSCAL, PSCAL, FCDNT, FxT 
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NY « FNY 

NZ » FNZ 

IF (NX.LT.l) GO TO 3Ci 
KPLOT * ABS(FPLOl) 

READ (iRFADfSOC) 

NN » 0 

11 NM - NM +i 

READ (1READ#51U F I T { NM ) » CGVO ( ) , P iu ( N?1 ) , ?<iO ( NM ) , P30 ( Nn ) f 

1 atTAO(NM)# STF 1P0(NM)>FHALF(NM) 

IF (FHALF(NM) .NE.O..ANO.NM.L1 .3) GO TG 11 

FHALF(3) - 0. I 

READ (IPtA0,500) 

read (IREA0>51C) FMACHy YA# AL , coo 
YAW = YA/RAD 

ALPHA ■ AL/PAD 

CALL GEOM (ND#NC»NP,ZS»XS#YS#XLE»YLt:f 5 lDP1>TRAIL> XP, YP, 

1 SWFtPl, SWtcPL>SwfcLP#DlHE01#DXHED2^DHED, 

2 XlEO#CHOKDw#ZTlP, iSYMO^KSYM) 

ISYM » ISYMO 

IF (ALPHA.NE.O. ) liYM » 0 
IF (KSYM.NE.O) YAW = 0. 

CYAW - CCS(YAW) 

SYAW « SIN(YAW) 

CA » CYAW+CQ$ (ALPHA) 

SA » C YAW-i-SINC ALPHA) 

IF (FCQNT.lt. 1.) GO TO vi 

READ (A) NX#NY,NZ,NM/M#K2#N IT 

MX • = NX +1 

MY • NY +2 

MZ * NZ +3 

DC 62 K»1,MZ 

READ (A) ( (G(l^ J>1 ),I*i>MX), J=1,MY) 

BUFFtR QUT(N3#1) ( G ( 1> 1 ) # G ( MX # M Y# 1 ) ) 

IF (LNIT(N3).GT.O.) GO TC i 

BUFFER QUT(N1,1) ( G ( 1, 1, 1 ) , G ( MX » M Y» 1 ) ) 

IF (UNIT(NI).GT.O.) GO TC 1 
62 CONTINUE 

READ (A) ( E0(K),K=K1,K2) 

REWIND N3 
REWIND NI 
REWIND A 

91 CALL CDQkC { N X> NY# NZ^ i<>i» YM » XI EC > ZT I P , XMA X r ZM AX, 

1 SY»iiCAL»i.CALZ#AX#AY»AZ# 

2 AC»A1#A2#A3#B0#B1#B2#B3#Z#C1»C2»C3) 

CALL SINGL ( NC # N Z # K S Y M# KT E 1# K TE 2# CHGR DC » 

1 SwEEP1#SWEEP2#SwEEP»0IHED1#0IhED2#D1HE0» 

2 Z1)#XLE,YLE#XC»X2»XZZ»YC»Y2#YZZ# 

3 Z#Cl#C2,t 3# El# E2# E3# EA# E5 # InO) 

CALL SURF (ND#NE#NC#NX#NZ# ISYM#KSYM# KTE1#nTE2, SCAL# 

1 YAW# AU# Z# Zj#XC# YC#SLQPT# TRAIL# XS#YS# NP, 

2 lTEl#ITE2#iV#S0#Z0#XP#YP#01#02#D3#X# Y# INU) 

IF (IND.EU.O) GC TO 291 

IF (FCQNT.GE.l. ) GU TC 101 
NM ■ 1 
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MT » 0 

CALL ESTIh 

IF (IQ.EQ.O) GO TO 1 
REWIND N3 
REWIND N1 

101 WRIT‘S (IWRIT^eOC) 

FCQNT = 0. 

MIT « FIT(NM) >NIT 

KIT » M I T 

IF (NM.GT.l, AND.FHALF(NM) .EU.G, ) MT ■ 10 
dIT • NIT 

KRES = (f'lT -NIT -a)/500 +E 

JPEi) » 0 

NR65) » 0 

COV » CDVC(NM) 

STRIP * STRiPC(NM) 

BETA ■ BETAO(NM) 


MX 

= NX 

♦ 1 

MY 

■ NY 

♦ 2 

MZ 

=« NZ 

+ 3 

KY 

■ NY 

+ 1 

K1 

« 2 


K2 

■ NZ 


IF 

(KSYM. cO.D ) 

GO TO iu3 

K1 

* 3 


K2 

• NZ 

+ 2 

103 LZ 

* NZ/2 +1 

IF 

(KSYM.Nt.O) 

IZ * 3 


WRITE (IWRITMOA) 

10^ FDRMAT(48H0INDlCATICiN uF LOCATION OF WING AND VORTEX SHEET# 

1 27h IN CCCRUlNAfE PLANE Y « 0./ 

2 27HG( ( IV(I#K)#K«K1#K2)#I»2#NX) ) 

DC 106 I«2#NX 

106 WRITE (IWPIT#65C» ( IV( I#K) »K=M#K2 ) 

WRITE (IWKIT#6GC) 

WRITE (IWRIT>112) 

112 FCRMATC49HCCH0RCWISE CtLL DISTRlbUTICN IN SQUARE ROOT PLANE# 

1 AMD MAPPEG SURFACE CODROINATES AT CENTER LINE AND TIP/ 

2 15H0 X #J5H ROOT PR0FILE#15H TIP PROFILE ) 

DO 11^ I»2#NX 

llA WRITE (IWRIT#610) AO ( I ) # SO ( I# L Z ) # SO ( I # KTE2 ) 

WRITE (lWPIT#ilt) 

lie F0PMAT(15H0 TE LOCATION »I5H POWER LAW ) 

WRITE (IWRIT#61C) XMAX#AX 
WRITE (IwRIT#60C) 

WRITE (IWRIT#11E) 

116 FDRiAT(46H0NQRMAL CELL DISTRIBUTION IN SQUARE ROOT PLANE/ 

1 15H0 Y ) 

DO 120 J«2#KY 

120 WRITE (IWRIT#610) BC(J) 

WRITE (IWPIT#122) 

122 FCRNATUSHO SCALE FACTDR#i5H POWER LAW ) 

WRITE nwRIT#tlC) SY#AY 
WRITE (IWRIT#60C) 
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WRITE (IWRIT#12A) 

12 A F 0 RMAT( 45 H 0 SPANWISE CtLL DISTRIBUTION AND SINGULAR LiSt/ 

1 15 HO Z , 15 H X SING , 15 H Y SING » 

2 15 h XZ , 15 H YZ X 2 Z » 

3 15 H YZ 2 ) 

DO 126 K»K 1 #K 2 

126 WRITE (iWRIT^eiC) 2 ( K ) # XC < K ) # YC ( K ) # XZ ( K ) » YZ ( K ) » XZZ ( K ) » YZZ ( K ) 
WRITE (iWRIT^lZb) 

128 FCRiATdSHO TIP LaCATiQNd 5 H POWER LAW ) 

WRITE (IWRIT# 610 ) ZMaX>A 2 
WRITE (IWRIT# 600 ) 

WRITE (IWPIT, 132 ) 


15H 


NZ 


132 FURNAT( 19H0ITERATIVE SGLUTiCN/ 

1 <r3HOSTRIP WIDTH EUR HORIZONTAL LINE RELAXATION) 

WRITE (lWRIT/610) STRIP 
WRITE (IWRIT,13^) 

13<> FORMATdSHO NX diih . NY 

WRITE dwRIT»6A0) NX,NY#NZ 
CALL SECOND(T) 

WRITE dWRIT^TOO) T 
WRITE dWRIT,13t) 

136 FORMATdSHO NACH NO *i6H YAW 

WRITE dwRIT>61C) FMACHfYA^AL 
WRITE (IWRITd36) 

138 FORNATdOHOITEFATiONdaH CORRECTION 

1 lt>H KL SIDUAl 

2 lOH CIRCULATNdOH REL FCT IdCH REL FCT 2,iOH REL FCT 3f 

3 ' lOH BETA dOH SONIC PTS) 


>li?H ANG OF ATTACK) 


dH 

^4H 


# 4 H 
»4H 


>4H 
f 4H 


K 

K 


141 NIT 
JIT 
PI 
P2 
P3 




= NIT +1 
• JIT +1 
«= PiO(NK) 

» P2C{NM) 

« P30(NH) 

IF (NlT.LE.iC) PI • 1, 

IF (NIT.LE.IO) P3 » 1. 

CALL HIXFLO 

IF dC.FC.O) GO TO i&l 
JO =0 

REWIND N1 
Rewind N2 
N * N 1 

Ni • N2 

N2 * N3 

N3 « N 

WRITE (lWklT>660) NIT>DG#IG»jG^K&>FR#IP#JR>KR>EU(LZ)# 


;UC®W'^^„'-’pO0R 

y^GE ® 


JRES 


Pl#P2#P3dFTAi NS 
+ 1 


- JRES 

IF ( JRES.EC.KRES) JRES = 1 
IF (JRES.NE.I) GO TO 143 
NR6S ■ NRFS 

CuUNT(NRES) = NIT 
RES(NRES) « FR 

143 IF (JIT. EC. KIT) 60 TO 251 

IF (MT.LT.KIT.AND.AES(OG).GT.COV.ANO.AaS(DG),LT.iO.) 


♦ 1 
-1 


GO TO 141 
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GO TO 161 

161 IF (JO.EQ.l) GG TO 1 
REWIND N1 
REWIND N2 
JO -1 

N » N3 

N3 - N2 

N2 * N1 

N1 « N 

GO TO 1^1 

161 Rate - 0, 

IF (NRES.GT.l) RATE = ( ABS (RES( NRES ) /RES(l) ) ) 

1 ♦♦U./(CCUNT(NRES ) -CQUNTd))) 

WRITE (IWRIT,162) 

162 F0RMAT(16H0 NAX RESIDAL 1,15H MAX RtSIDAL 2d5H WORK / 

1 15H REDOCTN/CYCLE) 

WRITE (IWRIT,67C) R E 5 U ) ^ Rc 5 ( NR t S ) ^ C OUN T ( NR ES )# RATE 
CALL SECOND(T) 

WRITE (IWRITWOO) T 
WRITE (IWRIT»6C0) 

DO 164 L»l#3 

BUFFER IN (Nl/1) ( 6 (1 L > # G ( MX, MY> L ) ) 

IF (UNIT(M).GT.O. ) 60 TO 161 
164 CONTINUE 

LX » NX/2 -H 

K » 2 

171 K * K +i 

IF (K.EQ.MZ) GO TO 191 
DC 172 J-1,MY 
DO 172 I'1,MX 
6(1, J,l) » G(I, J,2) 

172 G(I,J,2) = G(1,J,3) 

BUFFER IN (Nl,l) ( 6 ( 1, 1, 3 ) , G ( MX , MY, 3 ) ) 

IF (UNIT(Nl).GT.O.) GO TO 151 
IF (K.LT.KTEl.OR .K.GT.KTE2) GO TO 171 

11 » ITEKK) 

12 « ITE2(K) 

CALL VELO (K, 2, SV,SM,CP,X, Y) 

CHORD(K) - X(I1) -X(LX) 

CALL FDPCF ( 1 1, 12, X, Y , C P, AL, C HDRD ( K ) , XC ( K ) , 5C L ( K. ) , SCD ( K ) , SCM ( K) ) 
IF (rPLQT.GT. l.AND.K.GT.KTEl) GO TO 185 
WRITE (1WR1T,6CC) 

WRITE (IWRIT,102) 

182 FORMAT! 24H0SECT10N CHARACTERISTICS/ 

1 15H0 MACH NO ,15H YAw ,15H ANG OF ATTACK) 

WRITE (IWRIT,610) FMACH,YA,AL 
WRITE (IWRIT,184) 

184 F0RMAT(15H0 SPAN STATION, 15H CL , 15H CD , 

1 15H CM ) 

185 WRITE (IWRIT,61C) Z ( K ) , SC L ( K ) , SC D ( K ) , SC M ( K ) 

IF (KPLOT.LE.l) CALL CPLOT ( 1 1, 1 2, FM ACH, X, Y, C P ) 

IF (KPL0T.LT.1.0R.KPLQT.GT.2) GO TO l71 

CALL GRAPH ( I PLC T, I i, 1 2, X, Y, C P, T ITLE , FMACH, YA, AL, 

1 Z (K),SCL(K),SCD(K ),CHORDO, XSCAL,PSCAL) 
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GC TO 171 

191 CALL TQTFDR(KTEl,KTfci:^ CHGRO» JjCL^SCD# iCH, XC^ 

1 CL#CC1»C«P, CMR^Cr.Y) 

CDl • CYAW^CDl 

CD « CDO +C01 

VLDl « 0* 

IF (ABS(CDi).GT.l,L-fc) \/LOi = CL/COi 
VLD « 0. 

. IF (ABS(CL).GT. 1.6-6) vLD « CL/CD 
WRITE <I«RIT»600) 

WRITE (IWRIT>192) 

192 FORMATCaiHOWING CHARACTERISTICS/ 

1 15H0 ^ACH FO ,li>H YAW »16H AEG LF ATTACK) 

WRITE (IwRIT#610) FMaCH^YA^AL 
WRITE (IwRIT»19^) 

19A FGRMATdSHO CL dSH CD FORM #15H CD FRlCTiON , 

1 15H CD fl5H L/D FORM , 15H L/D ) 

WRITE (IWPIT,61L) C L^ CDl, CDO» CD# VLDl, VLD 
WRITE (IWRIT#196) 

196 F0RMAT(15H0 CM PITCH #15H CM ROLL #15H CM YAw ) 

WRITE {IWRIT#61C) CMP#CMR#CMY 
REWIND N1 

IF (KPLOT.lt. 1) GO TO £01 

CALL RPLOl ( IP LOT# MR ES# RES# COUNT# TITLE #FMACH# Y A# A L# NX# NY# NZ ) 

CALL THREEO(lPLC r#SV#SM#CP# X# Y# TITLE# YA# AL# 

1 VLD#CL#CD#CHORDG# XSCAL# PSCAL) 

IF (ID.EC.O) GO TO 151 
201 IF (ISTOP.EO.l) GO TO 301 

IF (FHALF(NM) .EC.O. ) GO TO 1 


NX 

■ NX 

♦ NX 

NY 

■ NY 

♦ NY 

NZ 

■ N Z 

♦ NZ 


CALL COORD ( N X# N Y# NZ, KS YM# XT EO# 7. T I P# XMA X# Z MAX# 

1 SY#SCAL#SCALZ#AX# AY#AZ# 

2 A0#A1#A2# A3#60#61#B2#b3# Z#C1#C2#C3) 

CALL SINGL ( NC# NZ # KS YM# K1 El# KTt 2# CHOkDO# 

1 SWEEPl#SWEEP2#iwEEP#DIHECi>0IHE0?#DIHED# 

2 ZS#XL£# YLE#XC#XZ#XZZ# YC# YZ#YZZ# 

3 2#Cl#C2#C3#£l#E2#Ei#£k#E5#lN0) 

CALL SURF (ND#Nt#NC#NX#NZ#ISYM#KSYM#KTtl#KTtf!#SCAL# 

1 YAW#Ao#Z#ZS»XC#YC#SLOPT#TRAlL#XS#\'S#NP# 

2 ITEl,ITE2#iV#S0#Zu#XP#YP#Dl#D2#03; X#Y# IND) 
IF (IND.EU.O) GO TO t9i 

CALL REFIN 

IF (lO.EO.O) GO TO 221 
REWIND N1 
REWIND N2 

NSMOC - -FHALF(NM) 

IF (NSMDO.lt. i) GO TO 211 
DO 202 N-1#NSMQ0 
CALL SMGO 

IF (lO.EQ.C) GO TO 221 
REWIND N1 
2t2 REWIND N2 
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211 

N 

6 

M 


N1 

S 

N2 


N2 

= 

N3 


N3 

■ 

N 


NM 

S 

NM 


MT 

s 

0 


GO TO 101 



221 

NX 

= 

NX/2 


NY 

K 

NY/2 


NZ 

B 

NZ/2 


CALL COORD (NX,Nr^ N2#KiYM» XTf ZTiPf XMAX# 2MAX, 

1 SYfSCAL# SCALZ# AX, AY, AZ, 

2 AO, A1,AZ,a3,6J,B1,P2,B3,Z,C1,C2,C3) 

CALL SINGl (NC,KZ,KSY'1,KTEl,KTfc2,CH’JK00, 

1 SwEtPi, SWLcP2# SWfcfcP,DIHEDi,DIHtD2,DIHcD, 

2 ZS,XLE,YLE,XC,XZ,XZZ, YC,YZ,YZZ, 

3 Z , C 1, C 2, C 3, Ei , b 2, L3 , E^, E t , I NO ) 

CALL SURF {ND,NE,NC,NX,nZ,ISYM,KSYM,K1E1,KTE2,SCAL, 

1 YAW,AO,Z,Z:>,aC, YC, SLOP! » TRAIL, XS,YS,NP, 

2 ITEl,rTE2,IV,SC,ZC,XP,YP,Dl,D2,D3,X,Y, IND) 

IF (IND.FC.O) GC TO 291 

GO TO 151 

251 K1 » KTEl ~i 

K2 - KTE2 +ITE2(KTE2) -NX/2 

DO 252 M*l,3 

WRITE (A) NX,NY,NZ,NM,Ki,K2,MT 

LD ZtZ K«1,RZ 

BUFFER IN (M,l) ( G ( i., 1, i ) , G ( MX , «Y, 1 ) ) 

IF (UNIT(Nl) .6T.0. ) GO TO 281 
262 WRITE (4) < (G(I, J,l), I»1,MX), J«1,MY) 

REWIND N1 

WRITE (4) ( FO(K ),K«R1,K2 ) 

ENDFILE 4 

252 CONTINUE 
REWIND 4 

CALL SSWTCHa,lSTQP) 

IF (ISTOP.EC.i) GO 7C 161 
JIT - 0 

IF (MT.LT.MIT.AN0.Af}S(DG).GT.CUV.AN0,A6S(DG).LT.10.) GO TO 141 
GO TO 161 
281 REWIND 4 
GO TO 151 

291 WRITE (IWR1T,600) 

WRITE (IWPIT,292) 

292 F0RMAT(24h0BA0 DATA, SPLINE FAILURE) 

GO TO 1 

301 IF (KPLOT.GT.O) CALL PLQ T ( 0. , 0. , 999 ) 

STOP 

50C FQRMAT(IX) 

510 FURMATC BFlO.6) 

530 F0RMAT(2OA4) 

600 FORMATUHl) 

610 FQRMAT(F12.4,7F15.4) 

620 FDRMAT( 8E15,5) 
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I 


630 FORMAT(1HO#20A4) 

6^0 F0«r1AF( le, 7115) 

65C FDRMAT<1X>32IA) 

660 FORHAK 1 10, El 5 . 5 , 3IA, fc 15 . 5, 3 1 5F10 . 5» in ) 

670 F0R'1AT(2E15.4,2F15.^) 

700 F0R»1AT115H0CQMPUT1NG TIME , FIG . 3 , lOH SECONDS) 
END 


C 


SUBROUTINE GEOh ( NO, nC, NP, Z S, X S, YS , XL E , YLE, SLOPT , TR A I L, X P, Y F, 

1 SWtEPl, SWEEP2, SWEE P,0IHtDl,DlHE02, OlHtU, 

2 XTEO>CHOkLO,ZTIP, lSYiO,KSYM) 

GEOMETRIC DEFINITION OF wING 

DIMENSION XS(ND,1),YS(N0,1),Z5(1),XLE(1),YLE( 1), 

1 SLOPT(l), rRAIL(i),XP(l),YP{i),NP( 1) 

IkEAO » 5 

IWRIT - 6 

KAO - 57.29577*^5130623 

READ (IRfcAD,5CC) 

READ ( IRE AD, 510) Z S YM , FNC , Sk LE P 1, S WE EP 2, S k E £ P , D IH E Cl , D iHtJ 2, t IHE D 

IF (FNC.LT.3.) RETURN 

KSYM » ZSYM 

NC « FNC 

WRITE (IWRIT, 2) 

2 FORMATdShO SwEEP(i) ,15H SWEEP(2) ,15H FINAL SWEEP , 

1 ■ 15H llHbD(l) ,15h DIHEu(2) ,15H FINAL DIHbD ) 

WRITE (IWRIT, 610 XL, YL, CHORE, THICK, AL 

WRITt (IWRIT, 610) $WEEPl,SWEf P2 , Sweep, EHED1,UIHEU2, DINED 


SWEEPl 

« SWEEPI/KAD 

SWE6P2 

« SWEEP2/RAD 

SWEEP 

■ SWEEP/RAD 

DIH6D1 

- DIHEDl/RAO 

DIHt02 

- 0IHtU2/RAD 

DINED 

* OIHED/RAD 

ISYMC 

» 1 

XTEO 

■ 0 • 

CHORDC 

« 0. 


K = I 

11 READ (IREAD,500) 

READ (IREAD,510 Z S ( K ) , XL, YL , C HQRD, TH I CK, AL , FS EC 
ALPHA « AL/RAD 

IF (K.GT.l.AN0.F:>tC.EO,O. ) GL TO 31 
( IREAO, SCO 

(IREAD,510 YSYM,FNU,FNL 
• FNU 
- FNL 

■ NU +NL -I 
( IREAD, SCO) 

(IREAD,510 TPL,SLT,XSING,YSING 
( IREAO, 500) 

I«NL,N 

(IREAD, 510 XP(1),YP(I) 


READ 

READ 

NU 

NL 

N 

kEAO 
READ 
READ 
DO 12 
12 READ 
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( 


L ■ iSL +1 

IF (YSYM.GT.O.) &u TO ii; 



READ 

( IREAD# SCO 


DO lA 

I-1,NL 


READ 

(IRFAD*510) VAL,DUrt 


J 

• L “I 


XP( J) 

■ VAL 

14 

YP( J) 

- DUN 


GO TO 

21 

15 

J 

= L 


DO 16 

1 «NL # N 


J 

= J “1 


XPU) 

- xp ( i ) 

16 

YPU) 

= “YP ( I ) 

21 

WRITE 

{ IWRIT/ 600) 


WRITE 

(IwPlT»22) ZS{R) 


22 FCRMAT(16H0PR0FILt A1 Z ■ >FiC.!?/ 

1 i^HO TE angle »l5h TE SLOPE 

2 l&H Y SING ) 

WRITE (IWkIT/610) TRL#SLT»XS1NG» YSl.MG 
WRin {IWRIT»2^) 

24 FORMAKl.^HO X f ISH Y 

00 2e 1»1#N 


p 15H 


) 


26 WRITE 
3i SCALE 
XLE{K ) 
YLfc (K ) 
XX 
YY 
CA 
SA 


(1WR1T>6Il) XP(I)>YPa) 

= ChORO/( XE (1) -XP(NL) ) 

• XL +(XS1n 6 -XP(NL) )»THlCK+SCALe 
= YL +(YS1KG -YP(NL ) )+THICK*SCALE 
» XP(M) -KXSiNO -KP(NL) »*TH1CK 
= YP(NL) +{YSiNG -YP(NL) ) ♦THICK 
« CQS(ALPHA) 

= SIN(ALPHA) 


32 


DC 32 I- 
XS( I,K ) 
YS(I#K ) 
SLOP! (K ) 
TRAIL (K ) 
NP(K) 
XTEO 
CHQRDO 


l^N 


* 5CALE+((Xpm -XX)4CA +TH iCK ♦ ( YP ( I ) 
- SCALt + (ThiCK + ( YP{ I) -YY)+CA -(XPd) 

* THICK+SLT -TAN( ALPHA) 

■ ThICK+TRL/RAL 
= N 

■ AMAXK XTEO^ XS( 1#K ) ) 

« AMAXKCHLROO# CHORD) 

IF (YSYH.LE.O..CR.ALFHA.Nt:.C. ) ISYHO « 0 
WRITE (IwPIT^52) ZS(K) 


52 FORMAT( 27H0S£CTI0n DEFINITION AT Z = #F10.5/ 

1 15hO XLE YLE ,15H 

2 15HTHICKNESS RATIOdbH ALPHA ) 

WRITE <iWRlT,61G) XL ^ YLj CHCRDdHlCK, AL 

K = K +i 

IF (K.LE.NC) GO TO 11 

ZO » .5A(ZS(1) +ZS(NO) 

IF (KSYM.NE.O) ZO « ZS(i) 

DO 62 K=1#NC 

62 ZS(K) - ZS(K) “ZO 

ZTIP « Z5(NC) 

RETURN 

500 FQRHAT(IX) 


X SING 


-YY)*SA) 

-XX)^SA) 


CHORD 
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510 F0(^MAT(6F10.6) 

600 FQRMAT(lHl) 

610 FORMAKFIZ.^, 7F15.4) 
END 


SUBROUTINE COCFD ( NX # N Y# NZ j K S Y K # XT bO , Z T I P / XMA X » ZM A X# 

1 SYf SCAL>SCALZ#AX# AY#AZ, 

2 AOfAl;AZ^A3;bO>Bl>bZ#B3«Z>Cl>CZ^C3) 

C SETS UP STRETCHIO PAKA30LIC AND SPANwlSc COORDINATES 

DIMENSION A0(l)>Al(l),A2(l)»A3(l)#b0(i)#bl( 1)^62 (l)f63(i)» 
1 Z(l),Cl(l),C2(i)#C3(l) 


DX 

m 

2. /NX 

DY 

m 

i./N Y 

KY 

X 

NY +1 

DZ 

B 

2./NZ 

ZO 

B 

1. -DZ 

K1 

B 

7 

K2 

s 

HI 

IF (KSYM 

.EQ 

.0) GO TO 1 

DZ 

B 

l./NZ 

ZL 

O 

0 . 

K1 

= 

3 

K2 

B 

NZ +2 

AX 

B 

ir 

• V 

AY 

B 

. 5 

AZ 

B 

. 5 

BX 

B 

0. 

BZ 

S 

0. 

XMAX 

B 

• 625 

ZMAX 

S 

.625 

SY 

B 

. 5 

SCAL 

K 

XTEC/ ( .500ul+ XMAX +X MAX) 

SCALZ 

m 

ZTIP/ U.OOOOOl+ZRAX ) 

V2 

=£ 

(0X/DY)*+2 

W1 

m 

SCAL /SCALZ 

k2 

X 

( Wl*0X/DZ)++2 

S73 

s 

SORT (73.) 

BBX 

= 

-BX+$0»T{3.t(7. + S73))/((l. + S 73 ) ♦ XM A X+*3 ) 

ABX 

a 

1. - BBX*SQRT( (7. + S73 ) / 12 . ) + XM AX + * 3 

CBX 

e 

(19. + S73 ) ♦XMAX+XMAX/12. 

ABBX 

B 

ABX ♦ BBX+(3.tC«X - A,*XMAX+XMAX)*XMAX*XMAX/ 
SORKCaX - XMAX^XNAX) 

DO 12 I* 

2»NX 

OD 

X 

(I -D+DX -1. 

B 

= 

1. 

IF (AeS(OD) 

•GT.XMAX) GO TO 13 

A 

X 

CBX - DD+DO 

AS 

= 

SQRT( A) 

C 

B 

ABX^AS + d3X+(3,YC8X - A . *0 D*CD ) +00 ♦ OD 

DO 

B 

ABX+DD + 88X*AS»DD+»3 

D1 

B 

AS/C 


60 



k 


I 


p 



D2 

m 

8PX-* + 19. + dD»uO) - 12.*00^*4)+D0/(A*C) 


GO TO 

14 


13 

IF (DD 

.LT.O. 

) 6 - - 1 . 


A 

= 

1. -{(OC -B + XMAX)/(1. -XF.AX))* + 2 


C 

■ 

A»%AX 


D 

S 

(AX +AX -l.)*(l. -A) 


DO 

X 

3+XFAX ♦ AFBX*(0r - P^XiAX)/C 


Di 

= 

A*C/((i. + 0)+AB6X) 


D2 

■ 

-(AX +AX)*(DD -6>^Xf1AX) 


1 


»(3. '►U)/((l. +D) + A + (1. -XKAX)* + 2) 

14 

AO(I ) 

a 

00 


Aid) 

a . 

. J)4C 1 /ox 


A2( I ) 

a 

01 + Dl 

12 

A3( I) 

a 

. i> + DX»02 


DC 22 

J-2»KY 


DD 

a 

(KY -J)*OY 


A 

a 

1. -DD»00 


C 

a 

A + + A Y 


D 

a 

(AY +AY -1,)*(1. -A) 


Di 

a 

A+C/((i. +D)+SY) 


t}0( J ) 

a 

SY^DO/C 


Bl( J) 

a 

.!?»D1/0Y 


b2 < J ) 

a 

01+DI+V2 

22 

B3( J ) 

s 

-AY»00*UY»( 3. +D)/((1. +D)i‘A) 


BBZ 

a 

-E2+S0BT( 3.»(7. + S73))/((l. + S 73 ) ♦ ZN A X+»3 ) 


A6Z 

a 

1. - BbZ+SURT((7. + S7a) /Ic. )+ZMAX^+3 


C6Z 

a 

(19. ♦ S73)*ZMAX*ZyAX/12. 


AB6Z 

= 

ABZ + 33Z+(3.»CBZ - 4.*ZrAX*Z^1AX) ♦ZhAX + ZMAX/ 


1 


SQRKCbZ - ZNAX + ZNAX) 


Do 32 

K-2#K2 


DD 

a 

(K -KD+OZ -ZO 


8 

a 

1. 


IF (ABS(DD). 

GT.ZNAX) GO TO 33 


A 

a 

CbZ - DD^OD 


AS 

a 

SCRT( A) 


C 

8 

A6Z+AS + 6BZ+(3.+CBZ - 4 . +0 D+D D ) +DD ♦ OD 


DO 

a 

AdZ^OC ♦ B8Z»AS»DD+43 


L 1 

a 

AS/C 


D2 

a 

BWZ+(CBZ*(-6.*CBZ + 19.#DD+D0) - 12 . *0 ) *0 J / ( A+C ) 


GO TO 

34 


33 

IF (DD 

•LT.O. 

) B ■ -i. 


A 

= 

1. -((DD -B*ZMAX)/(1. -ZMAX))i*2 


C 

a 

A + + AZ 


D 

« 

(AZ +AZ -l.)+(l. -A) 


DO 

a 

B + ZFAX ♦ AtJbZ*(00 - H*ZrtAX)/C 


DI 

a 

A^C/( (1. -R D) ♦AfleZ) 


D2 

a 

-(AZ +AZ)«(uD -b + Z(iAX) 


1 


♦(3. +D)/((1. +D)+A*(1. -ZHAX)*^2) 

34 

Z(K) 

a 

SCALZ^DO 


Ci(K) 

a 

. 5 + L l + Wl/ OZ 


C2(K) 

a 

D1+D1*W2 

32 

C3(K) 

RETURN 

EnO 

a 

. 5^DZ+02 
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subroutine SINGL ( nC » NZ# K Tfcl> KT £2 , CHORDO# 

1 SiKtfctPi#S*i(EfcP2»Sw£fcP^DIHfc:DijOlHcDZ#DIHtGj 

2 ZS>XLE,YLE,XC,XZ^XZZ>YC,YZ,yZZ, 

3 Z»Ci^C2#C3#El#h2#ta#t^#t5,IND) 

GENERATES SINGULAR LINE PQR SCUARE ROOT Tft ANSFQH Ma TlCjN 
DIMENSION ZS(l)»XLL(l)>YLE(l)>XC(i)#XZ(l)#XZZ(l)# 

1 YC(l)>Y2(l)»YZ2(l)»Z(i)»Cl(l)#C2(i)#C3(l)f 

2 El(l)»E2(l)it3(l)fEA(i)>E5(l) 

DO 2 K«1,NC 


EA(K) 

■ Q. 

2 E5(K) 

» 0. 

K1 

■ 2 

K2 

* NZ 


IF (KSYM.EO.O) 60 TO 11 
K1 -3 

K2 - NZ +2 

KTEi » 3 

11 DO 12 K»K1>K2 

IF (Z(K) .LT.ZSa ) > Kiel ■ K +1 
IF (Z(K) .LE.ZS(NC ) ) KTE2 » K 

12 CONTINUE 

B ■ CHOROO 

51 » TAN(SWctPi) 

52 » TAN(SwEEP2) 

T1 » TAN(OIHtCl) 

T2 » TAN(DIHED2) 

CALL SPLIE (1#NC,2S^ XLE#Ei#E2iF3>i»Sl,l,S2»0#0.. IND ) 
CALL INTPL {KTEl>KT£2/Z^XC#l>NC#ZS»XLE,cl#E2,63>0) 

Call intpl (ktei#kTc2>z»xz>i>nc#zs>c1#E2»E3>ea»o) 

CALL INTPL (,K.T£1»KTE2#Z#XZZ»1>NC^ ZS» E2» £3# lA» E5>0) 
call SPLJF (l,NC»ZS»YL£>El#E2#E3,l^Tl#l#T2j0#0.f IND) 
CALL INTPL {KTEl,KTE2^Z#YC#l^NC,ZSjYLEf El>E2,fc3#0) 
CALL INTPL (KTE1,KT£2,Z#YZ»1,NC,2S^ 61#E2»E3^cA,0) 
CALL INTPL (KTEl,KTE2#Z#YZZ,l>NC#ZS»f 2,t3,EA,E3,C) 

S » B+TAN(SwEEP) 

Si « 8+Sl 



S2 

X 

3 + S2 



T 

X 

8+TAN(DIhED) 



T1 

X 

B*T1 



T2 

X 

6 + T2 



XCt2) 

X 

a.’KXCIS) -XC(4)) +XC(i>) 



YC(2) 

X 

3.*(YC(3) -YC(A)) +YC(5) 



IF (KSYM 

.NE 

.0) GO TO 31 



N 

X 

KTEI -1 



DO 22 K* 

Kl# 

N 



11 

X 

(Z(K) -Z(KTE1))/B 



A 

X 

EXP( ZZ) 



XC(K) 

X 

XC(KTEl) +S*ZZ -(SI 

-A) 


YC(K) 

= 

YC(KTEl) +T + ZZ -(Tl -D + d. 

-A) 


XZ(K> 

X 

(S +(S1 -S)+A)/b 



YZ(K) 

X 

(T +(T1 -T)+A)/B 



XZZ(K ) 

X 

(SI -S)*A/(6+6) 


22 

YZZ(K) 

X 

(Tl -T)tA/(b»6) 


31 

N 

X 

KTE2 +1 
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piftlQltiAL PAGE 5V. 


i 


DL 32 K»N,K2 

ZZ » (Z(K) -Z (K,rt:2) ) /6 

A - EXP(-ZZ) 

XC(K) = XC(KTE2) +(S2 -A) 

YC(K) » YC(KTC2) ♦I^ZZ 4(T2 -T)*(l. -A) 

XZ(K) « 4(S2 -S)*A)/E 

Y2(K) » (T +(T2 -T>+A)/B 

XZZ(K) « -(52 -S)4a/(B+B) 

32 YZZ(K) ■ -(12 -r)4A/(34B) 

PETUSN 

t.ND 


C 

C 


SLBRPUTiNt SURF ( NO » N£, NC> NX# NZ# I S YM, K S YM » KTE 1 T B2# SC AL# 

YAW#A0,Z,ZS» XC# YC#SLUPT#TRAIL#XS# YS#NP# 

: nfcl#ITE2# IV#SCfZ0#XP#Ye#01#D2#D3#X#Y,lN0) 

INTEKPQLATES MAPPED WING SURFACE AT MESH POINTS 
INTERPOLATION IS LiNcAR IN PHYSICAL PLANE 


:NSION 


PI 

TYAW 

Si 

DX 

LX 

MX * 

MZ 

IVO 

IVl 

Du 2 K«1»M2 

ITEUK) 

ITE2(K) 

DO 2 I»1#MX 

IV(I#K) 

SO(I#K) 

K « 


SO(NE»l)#XS(ND#l)#YS(ND#i)#ZS(l)#SLOPT(l)#TRAiL(l)# 
XC(l)fYC( l)#AO(l) #Z(l)#ZO(n#X(l)#Y( D# 

XP( 1 ),YP( 1 ) »D 1(15 #L 2 ( 1 )#D 3 ( 1 )» 
IV(NE#l)»NP(l)#IThl(l)#ITE 2 (l) 

3 . 1 Ai 5926 S 35 b 97 y 
T An(YAk) 

• S+SCAl 
2 . /NX 
NX /2 +1 
NX +1 
NZ 43 

1 -ISYN -ISYM -ISYM 
-1 -ISYM 


IV(I#K) ■ -2 

SO( I#K) • 0, 

K » KThi 

K2 -1 

K2 ■ K2 41 

K1 » K2 -1 

R 2 « 1 . 

IF (ZS(K2) -Z(N)) 21#2b#23 

R2 » (Z(K) -Z5(Ki) )/( ZS(K2) -ZS(Kl)) 

R1 * i. -R2 

C • R14XS(1»M) 4R24XS(1#K2) 

CC » SQRT((C 40/SCAL) 

00 32 I-2/NX 

IF ((AO(I) 4.N4DX) .LT.-CC) 11 « I 4 

IF ((AO(I) -. 54 DX) .LT.Ct) 12 » I 

CONTINUE 
ITEl(K) « II 
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n E2(K) 
cc 

20(K) 

KK 

P 

N 

Q 

DO 42 I 
42 X(I) 
ANGL 
U 
V 


« 12 

= AO(12)/CC 

- Z<K) -7YAw*(XC(K) 

» M 

■ R1 

* NPCKK ) 

- SC)RT(XS(1^KK)/C)/CC 
»2#NX 

■ Q+AO( I ) 

= Pi +P1 

■ 1 • 

» 0 . 


+Sl^AO( I2)+A0(I2)) 


DO 44 I»1^N 

R * S0R1(XS(I^KK)+*2 + Y S ( i » »*. K ) • ♦ 2 ) 

IF (R.EG.C.) GD TO 4t 

ANGL » ANGL +A 1 AN2 ( ( U»Y S I I , KK ) - V ♦> S ( I # KK ) ) » 

1 (U*Xi>(I/KK) 4V + Yb ( i^KK ) ) ) 


U 

V 
R 

XP( I) 

YP ( I) 

GO TO 44 
45 ANGL 
U 

V 

XP { I ) 

YP ( I ) 

44 CONTINUE 
ANGL 
ANGLl 
ANGLE 
ANGLl 
ANGLE 
T1 
T2 

CALL SPLIF 


XS(1>KK) 

YS(l^KK) 

SQRT((R +R)/iCAL) 

R+CUS ( .S+ANGL ) 

R*S1N( .5 + ANGl ) 

PI 

- 1 . 

0 . 

0 . 

0. 

ATA^ ( SLUPTCKK ) ) 

ATAN( Y5 (l^KK) /XS(i# KK) ) 

ATAN( YS(N#KK) /XS(N>KK) ) 

ANGL -.5*(ANGLi -TRAIL(KK)) 
ANGL -.5 + (ANGL2 +rRAlL(KI<)) 
TAN(ANGlI) 

TAN( ANGL2) 

(l,N>XP,YP#Di,02»C3#i»Tl»i,T2,Cf0 


# InO) 


CALL INTPL (Il^I2,X#Y»i>^»XP^ YP.»D1#D2#C3#0) 

XI » .25>* XS< 1#KK) 

A » SLDPT(KK)+(XS(1>KK) -XI) 

B - 1,/(XS(1^KK) -XI) 

ANGL * PI +Pi 

U » 1. 

V * 0 • 

N - II -1 

DO 52 I=2^M 

XX ■ .5ASCAL+X(I)++2 

D » 3+(XX -XI) 

YY ■ YS(l^KK) ♦A+ALOG (D) /O 

R » SQRKXX++2 +YY**2) 

ANGL ■ ANGL +aTAN2((U*YY -V+XX)#(U+XX +V+YY)) 

U « XX 

V « YY 

R = SCRT { ( R +R) /SCAL ) 
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5 ? 


54 


62 


71 

72 


74 


76 


81 

82 


Y(I) = R + SIN( . ':+MNtiL ) 

A - SLOPrCKiO^lXSlSfl-K) -XI) 

a * i./(XS(N,KK) -XI) 

ANGL ■ 0. 

U » 1 • 

V ■ 0 • 

=12+1 
Du 54 I»M#NX 


XX « .5 + SCAL4X I)» + 2 

D - B+(XX -XX) 

YY * YS(N»KK) +A'*ALaG(D) /D 

R » SQPT(XX**2 +YY^+2) 

ANGL ' ANGL ♦ Al AN2 ( ( U+Y Y -V+XX),(U*XX 

U ■ XX 

V = YY 

P • SQPTUP +R)/SCaL) 

Y{I) » R*S1N( , 5 + ANGL ) 

0 » P*C+CC*CC 

00 62 I»2#NX 

S0(I»K) - S0(I»K) tQ+Yd) 

IF (KK.E0.K2) GC TO ?1 
KK » K2 

P » R2 


GO TC 41 


LO 72 I«ilX2 
IV ( I»K) » 2 

t*- » 11 -1 


♦ V + YY) ) 


00 74 l«2fM 

Z2 • 2(K) -TYAW*( XC (K ) +S 1 ♦ AO ( I ) ♦ AO ( I ) ) 

IF (ZZ.Gt.ZO(KTFl) ) IV(I>K) «= IVO 
CONTINUE 

rt -12+1 

DO It I*f",NX 

ZZ - Z(*<) -TYAW+(XC(K) +S1AA0(I)»A0(1) ) 

IF ( 2Z.GF.Z0( KTtl ) ) IV(i»K) = I VC 
CONTINUE 

K2 ■ K2 -1 

K » K +1 

IF (N.Lfc.KTEZ) GO TO 21 
K1 -2 

K2 ■ NZ 

IF (KSYM.tQ.O) GO TC 61 
K1 - 3 

K2 » NZ +2 

DO 82 I»2#NX 

ZZ « Z(K) -TYAW»(XC(K) +S1»A0( 1)AA0( I) ) 

IF (ZZ.LL.ZS(NC) .AND.ZZ.G£.ZC(NTF1) ) 1V(I,K) ■ IVO 

CONTINUE 

K » K -t-l 

If <N.LL.N2) GO TO ei 
N ■ K TE2 

IF (YAW.LE.O.) go TO 93 
10 ■ ITEl(KTt2) ♦! 

00 92 I»1C#LX 
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N * N +i 

92 ZOtN) = Z(KT62) -T ( XC ( Kl £2 ) + S AC ( I ) ( 1 ) ) 

93 I - ITEKKTtl) 

ZO(KTfcl-l) « Z(KTtl-l) -TrAWt{XC(KTEl-i) +S1+A0 ( I ) ♦AO( I ) ) 
Z0(N + 1) - Z(KU241) 

DO i02 K»Ki#K2 
DO IGA I«2,NX 

IF (IV(I,K).GT.C) GO TO lOA. 

IF av( I + l/K + i) .GT .0.0k. I\/( I-1»K + 1) .6T.0) IV(1»K) ■ IVl 
IF (IV(I+l>K-i).GT.O.QR.IV{I-l#K-i).GT,G) IV(I#K) = IVi 
lOA CONTINUE 

102 IF (SO(LX>K ) .LT.l.E-05) i\/(LX#K) « 0 
IF (KSYM.tQ.O) FtTURN 
DO 112 I-2^NX 

112 S0(I#2) » 3.*(S0(I^3) -SG(i»A)) +St(I>t)) 

RETURN 

END 


C 


SUSRCUTINE ESTIM 

INITIAL ESTIMATE OF KcOUCED POTENTIAL 


COMMON 


1 

2 

3 

A 

5 

6 
7 


G{193,26^A)^S0(193#35)/t0(131)#Z0(131)^ 

IV(193,3i3),IT£i(33)>ITE2(35), 

A0(193),Al(i9 3),A2U93)# A3( 193)^ 

B0(25),ei(25),82(26),83(26), 

Z{35)#Ci(35)>C2(35)»C3(35)# 

XC(35)>KZ(35)^XZZ(36)>YC(35)^YZ(3:')/YZZ(35)f 
NX>NY^NZ>KrEl»KTt 2f I SYM> RSY M » :iC ALf SC AL Z > 

YAWf CYAWf SYAw# ALPHA,CA, SA, FMACH# N1, N2» N3» ID 


12 


21 


MX 



m 

NX 

+ 1 

KY 



= 

NY 

+ A 

MY 



■ 

NY 

+ 2 

MZ 



s 

NZ 

+ 3 

DO 

12 

I- 

l, 193 


DO 

12 

J = 

ip 26 



DO 

12 

K- 

IpH 



G ( I« Jj 

.K) 

m 

0. 


K 



m 

1 


DO 

22 

I» 

ZpHX 




G( I 

#KY+1,1) 

■ 0. 


IF 

(IV(I^K) 

.LT.2) GO TO 

22 

DSI 

m 

S0(I+1#K) - 

SOU-i^ K) 

DSK 

X 

50{1>K+1) - 

S0( dK-i ) 

SX 

m 

Aid )*DSI 


SZ 

X 

C1(K)+DSK 


FH 

X 

AO(I )^AC( I) 

♦ S0( dK)+SO( dK) 

H 

3 

i./fh 


AZ 

X 

-A0( I )>»XZ{K) 

-SO(dK)^YZ(K ) 

BZ 

X 

-A0( I )»YZ(K) 

+SC(I>K)+X2{K) 

HZ 

X 

Az+sx -ez 

+FH+SZ 

FYY 

X 

1. fSXtSX 

+H+HZ+HZ 

FXY 

X 

SX +H+AZ+HZ 
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V » S A*^0( I ) -CA + i>0( I, K) 

U » CA*AC(l) +bA+SO(I#K) 

W • SYAW +CA+XZ(K) +SA»YZ(K) 

G(I>KY+1>1) = G(I#KY-1#1) 

i ♦(V + (l. -h»8Z»HZ) -U*FXY -w*HZ ) /( FYY + BUkY) ) 

22 COSTINUE 

bUFFtk 01T(N3>1) (0(1, i# i)>G( ^X,f*Y, i) ) 

If (UNIT(N3) .GT.O. ) GQ TO 41 

BLFFtft QUT(Nlfl) ( G ( 1 , 1, i ) , G ( hX , MY, 1 ) ) 

IF (UNIT(NI) .GT.O. ) 00 TD 41 
K « K +1 

IF (K.LE.MZ) GO TO 21 
K1 = KTEI -1 

K2 • KTF2 +nt2(KrE2) -NX/2 

UO 32 K«K1,K2 
32 eC{K) » 0. 

l(J * 1 

RcTJPN 

41 10 » C 

ktTlJPN 
lNO 


SU3R(]UTlNt MlXFLU 

C SOLUTION OF EQUATIONS FOR hIXbC SUBSONIC AND SUPERSONIC FLOW 

C USING ROTATED DIFFERENCE SCHEME 

common G(193,26,4),S0(143,3&),E0(131), Z0(131) , 

1 IV(193,35),ITE1(35),ITE2(35), 

A0(193),Ai(193),A2(193),A3(193), 

8C( 26), 81(26), B2(2fa),E3(26), 
Z(3t),Cl(36),C2(35),C3(35), 

XC(35),XZ(35),XZZ(36),YC(35),YZ(35),YZZ(35), 
NX,NY,NZ,KTei,KTE2,ISYM,KSYt1,SCAL,SCALZ, 
YAW,CYA*,SYAk,ALPHA,CA,SA,FMACH,N1,n 2, N3,IQ 
COMMON /FLO/ S TR IP , P 1, P 2, P 3, B E TA , FR, I R, J R,X R, D G, IG, JG,KG,NS 
COMMON /SkF/ Gia( 193, 26),GK2{193,26), 

SX(193),SZ(193),SXX(193),SXZU93),SZZ( 193), 
R0(193),F 1(193),C (193),0( 193), 

G10(26),G20(26),G3C(26), 64C (26),G1(26) ,G2(26), 
il,12,K,L,N0,LX,MX,KY,MY,Tl,A40,Ql, Q2, T YAW, SI 


LX 

= 

MX/2 +1 

MX 

■ 

NX +1 

KY 

S 

NY +1 

MY 

a 

NY +2 

TYAW 

a 

SYAW/CYAW 

SI 

a 

. 5+SCAL 

DX 

a 

2. /NX 

T1 

a 

OX + CX 

AAG 

3 

1./FMACH++2 

Q1 

a 

2. /PI 

02 

a 

1./P2 

FR 

a 

D • 
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IR - 0 

JR » 0 

KR » 0 

OG * 0. 

IG - 0 

JG =0 

KG » U 

NS * 0 

K1 ■ E 

IF (FMACH.&E.l. ) K1 - 3 
K2 - NZ 

IF (KSYM.EQ.O) GO 70 1 
K1 « 3 

K2 ■ NZ +2 

IF - AliS( .5 + S7RIP + NX) 

L - F 

IF (L.EQ.NX/2) L = L -i 

11 - LX -L 

12 « LX +L 

IF (L.EQ.O) 12 « LX -1 

00 2 L=1^3 

BUFFER IN (Nl,l) (G(i»lU )#G(nX^("Y,L) ) 

IF (UMT(Nl) .GT.O.) GO TO 101 
2 CONTINUE 
DO 4 J«l>hY 
DO 4 i=l#NX 
G(I>J»4) ■ G(Ii J»l) 

GKl(i»J) - G(I#J»1) 

4 GK2(1#JJ - G{I»J»1) 

K « 2 

L -2 

NO “ KTEl -1 

IF (H.EO.Kl) GO TO 21 

buffer 0UT(N2>i) ( G ( 1 , 1, 4 ) , G ( MX , MY, 4 ) ) 

IF {UNIT(N2) .GT.O. ) GO TO 101 
BUFFER IN (Nl,l) (G(1,1,4),G(MX,MY,4) ) 

IF (UNIT(NI) .GT.O.) GO TO 101 
IF (KSYM.EQ.O) GO TO 51 

1 - LX 

OSI = ;.0(I+1,3) -S0(i-1,3) 

DSK « son, 4) -S0(I,2) 

SX(1) - A1(1)*0SI 

SZ(I) - C1(3)+0SK 

R = AMINOd, 1V( I,K) ) 

J ■ KY 

DO 12 M»2,KY 

YF • BO(J) +S0(I>3) 

H » R/(l. -R +YP+YP) 

AZ • -YP + YZO) 

BZ » YP+XZ(3) 

A » H + A2 + AK1) 

B * (H+(8Z -AZ^SX(D) -SZ ( I ) ) tBl ( J ) 

DGI - G(I+1,J,3) -G(I-1,J,3) 

DGJ - G(I,J+1»3) -G(I,J-1>3) 
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G(I,v'.2) = 

GK1(1,J) - 

G(IfJ#l) « 
GK2(i#J) - 
P = 

12 J 
J 

G(I#Ji?) - 
6K I (I , J ) * 
G(I#J»1) - 

GK2(1,J) « 


n « , 

DO H = lfM 

I 

GO ro 16 

15 1 

16 DSI 
DSK 
SX(I) 

SZ(I) 

DO IP J=2>KY 

YP 

H 

AZ ■ 

ez 

s 

A * 

B 

IP 

IM 

DGI 

Co J ■ 

G(I,J,2) = 

1 

GKl(lfJ) - 
Gd^J^l) - 
18 GK2(I»J) = 

J • 

G(I,J,2) » 

1 


GK1(1,J) = 

IF (l.LT.LX) 
lA CONTINUE 


6(r»J»A) +(A>QG1 -B + DGJ) /CK3) 

G(d J,2) 

3.>MG(dJ>2) -G(I>J,3)) +G(dJ#A) 

G(I# 1) 

1 . 

J -1 
KY +1 

G(dJ>A) ♦(A*DGI -6*DGJI'/Cl(5) 

G( I» J# 2) 

3.+(G(I#d#2) -G(i>J,3)) +G(I»J»4) 

6( I» J,l) 

NX/2 -1 

LX -II 

LX ♦!! 

S0(l+1>3) -3G(i-i#3) 

SO(I»A) -S0(I>2) 

Al( I )»DSi 
C1(3 )^CSk 


HO(J) +S0(If3) 

1./(aO(I )+AO(I) +YptYP) 

-AC(I)*X^(3) -YP*YZ(3) 

-A'j(l)*YZ<3) +YP*XZ(3) 

dIGN(l.,AZ) 

H+A8S( AZ)tAl(I) 

(H+(BZ -AZ + 3X(D) -SZ( i) )+81( J) 

1 -HFIX(S) 

I -ir-IX(S) 

G(dJ#A) — G(IN»J>A) 

G(I>J4!#3) -G(I#J-1>3) 

(C1{3)*G(I»J#A) +A*(G(IP# J#2) +DGI) 

(CK3) +A) 

G( I, J>2) 

3.»(G(I#d»2) -G(I,J#3)) iGd/J^A) 

Gd, J,l) 

KY -d 

( Ci(3)*G( d J# A) +A+(GdP,J,2) +DGI) 
(CK3) +A) 

Gd, 0,2) 

GL TO 15 


-6+D6J) / 


-B+DGJ) / 


GU TO 51 

21 BUFFt*? QUT<N2,1) ( G d, i, A ) , G ( NX , M Y, A ) ) 

DO 22 0=1, MY 

GiO(J) - Gd2,0,2) 

G20(J) » Gd2-l,J#2) 

G30(0) - G(I1,0,2) 

22 GAO(J) = Gdl + 1,0>2) 

00 32 1=2, NX 

DSI « S0d + 1,K) -S0( 1-1, K) 

DSK » SOd,K-fi) -S0(l,K-i) 

DSII = SGd + l,K) -$Od,K) -S0d,K) +SC(i-l,K) 
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1 

DSKK 

m 

+ A3( I )'^DSi 

SOn,K+l) -S0(1#K) -SO(I#K) +SC(I#K-1) 

1 

OSIK 

m 

+C3(K)^DSK 

S0(I+1#K+1) -S0(I-1#K+1) -S0(i+1#K-1) +SO(I 

SX(I) 

m 

Aid )*DSi 

SZ(I) 

m 

CKK )+05K 

SXX( I) 

= 

A2d )*DSIi 

SZZ ( I ) 

m 

C2(K)*0SKK 

32 SXZ(I) 

s 

T1+A1(1)+C1(K)+DSIK 


IF (12.GT.il) C/JLL YSWfccP 
IF ( UNIT(N2) .GT.O. ) GO TU lOi 

IF (K.LT.K2) BUFFER IN (M,l) ( G ( 1» 1# ^ > G ( M X » hY , 4) ) 

IF (I1.GT.2) CALL XShttP 

IF (UNIT(Nl) .GT.O. ) GO TC 101 

IF (t< .NE.KTE2.CP . YAw.LE.C . ) GO TQ t)i 

10 • ITEl(K) ♦! 


DU <»2 I = I(;,LX 
M ■ NX +2 -I 


£ ■ G(h#KY,2) -G(i»KY#2) 

NO « NO ♦! 

42 EO(NO) = EC(NU) +P3^(E -EO(Nu)) 

51 IF (K.EQ.K2) GO TO 61 


DO 52 J»1>MY 
DO 52 I«l>MX 


G( I, J,l) 
G(l,si*Z) 
G(I, J#3) 
52 G(I,J,4) 


G( I>Jj 2) 
G ( I > J # 3 ) 
G(I# J.4) 
G ( I # J j 1 ) 


K * K -n 

GO TO 21 

61 DO 62 L=«2»3 

BUFFER QUT{N2»1) ( G (1, 1, L ) MX# NY, L ) ) 
IF (UNITtNE) .GT.O. ) GO TO 101 

62 CQNTINUF 

FR ■ 1.2+FP/AAO 

10 = 1 

RETURN 

101 10 - 0 

return 

END 


SUBROUTINE YSWEtP 
C ROW RELAXATION 

CCMHON G(193#26#4),S0(193#35)#EC(131), Z0(131) # 

1 i V(193#35)#iT£l{35) #ITE2(35), 

2 Ag(193), A1(193)#A2( 193)» A3(193)# 

3 B0(26)#bl(26)#B2(26)#B3(26)# 

4 Z(35)#C1(35)>C2(35)#C3{35), 

5 XC(B3)>XZ(35)#XZ2(35)#YC{35)» YZ(35)#YZZ (35), 

6 NX»NY#NZ#KTE1#KTE2#ISYM#KSYH#SCAL#SCALZ# 

7 YAw#CYAW#$yAw#ALP HA # C A# S A # F N aC H#Ni#N2» N3#I0 
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I 


f 


CDM.>1CN/FlO/ 

COMrtGN/SwP/ 


1 

2 

3 

<» 


STRlP>Pl>P2^P?#BbTA>FR>IR# JRfKR^OGf IG/ J6^KC->NS 
&K1(193,25),GK2(193#26)^ 

SX(193)»SZ(193)fbXX(i93),SXZ(i93)>SZ2(193)# 
RO(193)^Ria93)>C (193)^0(193)# 

G10(26)#G20(25)#G30(26)#G90(26)#G1(26)#G2(26)# 
Il#I2#K#L#N0#LX#MX#KY#Mr#Tl#AA0#Qi#Q2# T YAW# SI 


12 


31 


J1 

m 

2 


IF (FMACH.GE 

.1.) J1 = 3 


C (11- 

1) 

0. 


D( 11- 

1) 

0, 


DO 12 

I-Il# 12 


Rod) 

M 

1. 


Rl( I ) 

U 

1. 


GKK 1 

#1) 

» G(1#1»L) 


GKld 

# Jl-1 ) 

- G(I#Jl-l#L) 


J 

X 

J1 


13 

m 

12 


BC 

m 

-T1 + 6K J)+C1(K) 


DO 32 

I-Il# 13 


AB 

s 

-ll^AK 1)+Bl( J) 


AC 

s 

Tl^AKI )»C1(K) 


YP 

X 

S0(1#K) +60(J> 


A 

X 

i. -RG(i) *AC(I)+AO(i) 4YP+YP 


h 

- 

R 0 ( i ) / A 


FH 

X 

RO(I )*A 


P 

X 

A0( 1 ) ♦( A.+YP+YP -FH) 


0 

X 

YP^(4.»A0( 1)'*A0( I ) -FH) 


A 

= 

XZ(K)*XZ(K) -YZ(K)*YZ(K) 


B 

X 

(XZ(K) +XZ(K) )^YZ(K) 


AZ 

X 

-A0(1)+XZ(K) -YP+YZ(K) 


BZ 

X 

-A0(1)4YZ(K) +YP*XZ(K) 


cz 

X 

HWH+(P*A -043) -AC(I)4XZZ(X) 

-YP*YZZ (K) 

D2 

X 

H4H4(U4A +P49) -A0(I)4YZ2(K) 

+YP4XZZ (K) 

DGI 

X 

G(I+1#J#L) — 3(1— 1#J#L) 


DGJ 

X 

G(i#J+l#L) -GKi(i#J-l) 


DGK 

X 

G(I#J#L+1) -GK1(1#J) 


DGII 

X 

G(I+1#J#L) -G(I#J#L) -G(I#J#L) 

+A3(I)*DGi 

+G(1-1# J#L) 

DGJJ 

L 


6(I#J+i#L) -G(1#J#L) -G(1#J#L) 

-B3( J )40GJ 

♦ G(I# J-1#L) 

DGKK 

L 

= 

G(I#J#L+1) -G(I#J#L) -G(I#J#L) 

■►C3((0 4DGK 

+ G(I# J# L-1) 

DGIJ 

L 

X 

b(I+l#J+l#L) -G( I-l#J+i#L) 
-G(I+1# J-1#L) 4G(I-1# J-1#L ) 


DGIK 

L 

X 

G(I+1#J#L+1) -G( I+l# J#L-1) 

-G( I-l# J#L-H) +G(I-1# J# L-1) 


D6JK 

L 

X 

G(I#J+1#L+1) -G(1#B-1#L+1) 

-G(I# J+1#L-1) +G(I# J-1# L-1) 


GX 

X 

Al( 1 )4DG1 


GY 

X 

-fll( J ) 4DGJ 


U 

X 

GX -SX(1)+GY +CA4A0(I) +SA4YP 


V 

s 

GY +SA4A0(1) -CA+YP 


w 

X 

RC(I)4(C1(K)40GK -SZ(1)+GY +SYAW 
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,-t (M ro sT 


AU - U +W+AZ 

AV « V +W+BZ 

QXY » H^(U+U +V^V? 

QQ ■ GXY 

AA • DIM( AAO# .2+UQ) 

HZ • AZ+SX( n -BZ +FH*SZ( i) 

FXX * 1. +H>»AZ*AZ 

FYY - 1. +SX(1)+SX(I). +H*HZ*HZ 

FXY » SX(1) +H^AZ*HZ 

BV - AV -AU*SX(I) -FH+W*SZ(i) 

UU * H^Al+AU 

VV ■ H^BV+BV 

WW » FH + iK + W 

UV » H+AU+8V 

= AU + V. 

VW ■ BV+W 

AXX ■ RKl )♦( FXX + AA ~LU) 

AZZ » FH+AA -W'w 

AXZ = (RKI) +RK1) )♦( AZ + AA -UW) 

R --(AXX^SXXd) +A7Z*SZZ(I) + A X Z* S XZ U ) > +G Y 

+Tl+( AA*(CZ*GX +(UZ -SX(I)+CZ)+GY) 

-HA(CA*( AU*AU -AV*aV) +(SA +SA)^AU»AV 


AXT 

AYT 

AZT' 

A 

AXT 

AYT 

AZT 

IF (CQ.GE 

AXX 

AYY 

AZZ 

AXY 

AXZ 

AYZ 

BP 

BM 

B 

R 

1 

GO TC 35 
33 NS 
S 

IH 

IMM 

AXX 

AYY 

AZZ 

AXY 

AXZ 


-UXY+(U+AC ( I ) +V*YP 

+(W +W)*( AC(i)AAZ +YP+BZ))) 

(Ca+XZZ(K) +SA*YZZ(K)) -w+w^(U^CZ +V^DZ)) 

- A6S( AU^AK 1) ) 

■ ABS(BV^E1( J) ) 

= AbS(FHtWACi(K) ) 

■ RO(l)tBETA + AA/AMAXl(AXT, AYT^AZT^d. -RO(I))) 

» A + AXT 

■ A + AYT 
« A^AZT 

AA) GO TO 33 

■ AXX*A2(I) 

« (FYY+AA -VV)*BZ(J) 

- AZZ4CE(K) 

= -Rl( I)>MFXY^AA +UV) + (AB 4AB) 

• AXZ4AC 

« -P1(I)4(HZ+AA +VW)4(BC +BC) 

- AXX 
« AXX 

« -AXX -AXX -Qi*(AYY +AZZ) 

- AXX+DGII +AYY+D6JJ +AZZ*DGKK 
♦AXY4DGIJ ♦AYZ^CGJK +AXZ40G1K +R 

* NS +1 

- SlGN(l.#b) 

« I -IFIX(S) 

« IM -IFIX(S) 

= UU+A2(1) 

- VV+B2(J) 

» WW*C2(K) 

« ti.+StUV+AB 

= 8.+S4UW4AC 
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AYZ 

B 


exx 

■ 

IFXXfQQ -UU)*A2(I) 

BYY 

= 

(FYYtOC -V\/)+B2(J) 

B2Z 

m 

(FH+QO -WW)*C2(K) 

BXY 

3 

-(FXYtQQ +UV)*(AB +A6) 

BXZ 

m 

(AZ*00 -Uw)»(AC AAC) 

BYZ 


-(HZ+QQ +VW)t(ttC +80 

AQ 

M 

AA/QO 

DcLTAG 

= 

BXX^DGII +BYY+OGJJ +BZ2+DGKK 

1 


♦ 6XY + DG1 J +tiYZ + CGJK +3XZ*DGIK 

DGir 

= 

G(I,J,l) -G(iM>J,L) “GdPiJ^U +G(iMM,J,L) 

1 


+A3(I)+0G1 

DGJJ 

■ 

G(I,J,L) -G(I,J-1,L) -G(I,J-1,L) +GK1(I#J- 

1 


-B3( J)*DGj 

DGKK 

3 

G(I#J»L) -G(I,J,l.-l) -G(I,J,L-i) +GK2(I#J) 

1 


+C3(K)+0GK 

DGIJ 

= 

G(I/u#L) 

1 


-G(I»J-1>L) +G(Ih#J-l#L) 

DGIK 

B 

G(I^J/L) -G(I»J#L-i) 

1 


-G(If*#J,L) +G( IK, J,L“i) 

DGJK 

= 

G(T,J,L) “G ( I , J , L""l ) 

1 


-G(I,J-1,L) +G( I, ) 

GSS 

m 

AXX4DGII +AYY+DGJJ +AZZ+CGKK 

1 


+AXY+DGIJ +AYZ+DGJK +AXZ+DGIK 

B 

m 

.5+(AC -l.)»(AXX +AXX +AXY +AXZ) 

riP 

= 

AO+bXX -(i, -S)+B 

BM 

■ 

A0+6XX -(1. +S)+B 

B 

c 

-AC+(SXX +8XX +Q2*(BYY +8ZZ)) 

1 


♦(AC “1 . )+(2.+( AXX +AYY +AZZ) +AXY +AYZ 

P 

s 

(AO -l.)*GSii +AQ + DtLTAG +R 

35 IF (ARS(R) 

• 

LE.ABS(FP)) GO TO 37 

Fk 

3 

R 

IR 

B 

I 

JR 

3 

J 

KR 

S 

K 

37 k 

3 

R -AYT+(GR1( I, J*l) -G(I,J-1,D) 

1 


-AZT*(GK1(I, J) -G(I,J,L-D) 

B 

3 

e -AXT -AY7 -AZT 

BM 

= 

6M +AXT 

B 

B 

i./(8 -BM + Cd-in 

C ( I) 

= 

B+BP 

32 D<I) 

3 

3*(P -eH + Od-D) 

CG 

= 

c. 

I 

3 

13 

DCj 42 M-Il 

» 

13 

CG 

8 

0(1) -C(I)+CG 

IF (ABS(CG) 

,i.e,Aes( 06) ) GO TO 43 

DG 

= 

CG 

IG 

8 

I 

JG 

3 

J 

KG 

3 

K 

43 GK2(i#J) 

= 

GK1(I,J) 

GKK I#J ) 

8 

G(I, J,L) 

G( I» J,L) 

3 

G(I,d,L) -CG 


+AXZ) 
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i 


42 I =1-1 

J ■ J +1 

IF (J -KY) 31,S1»51 

51 IF (12.GT,ITC-2{K) ) 13 • lTc2(K) 

IF { ITE2(K) .EO.MX) 13 * LX 

DO 52 1-11,13 

LV - lABSa -IABS(IV( I,K) ) ) 

RO(I) - AHINO(LV, XAB3( IV{ l,h) ) ) 

52 RKI) = LV 

GO TC 31 

61 N -NO 

I = LX +1 

IF (H.LT.KTEl.OR.K.GT.KTfcZ) GO TO 71 
10 » NX +2 -13 

DO 62 1-10,13 

A -1. -RO(I) +AO(I)^ACHI) -»S'H 1,K»+S0( I,K) 

H - R 0 ( i ) 7 A 

FH = R0(I)*4 

AZ ■ -A0( I)^X2(K) -S0( 1,K)»YZ (K ) 

BZ « -A01I)+YZ((<) +S0(I,K)*XZ(K) 

HZ - AZ^SXU) -aZ ♦FH>*SZ{J) 

FYY * 1. +SX(1)+SX(I) +H+HZ+HZ 

FXY - SX(1) tH+AZ+HZ 

DGI « G(I + 1,KY,L) -G(I-1,KY,U 

DGK - G(I,KY,L-H) -GK2(1,KY) 

V = SA4Ao(I) -CA+SC(I,K) 

U ■ A1(1)+DG1 +CA+A0(]) +SA*SO(J,K) 

W « C1(K)+DGK +3YAW +CA+XZ(K) +3A+YZ(K) 

62 G(I,KY+1,L) « G(I,KY-i,L) 

1 +(V*(1. -H+BZ*HZ) -U*FXY -W»HZ ) /( FYY+8 1 ( K Y ) ) 

I = 10 

IF ( lO.NE.lTcKK ) ) GO TO 71 
E » G(I3,KY,L) -G(i0,KY,L) 

NO * NO +1 

EO(NO) = EO{NO) +P3+(E -EO(NO)) 

N -NO 

71 IF U.LE.Il) RETURN 
I • I -X 

E =0. 

IF ( 1V( I,K) .NE.l) GO TO 77 

ZZ - Z(K) -TYAW+(XC(K) +S1 + AO ( 1 ) +A0 ( I ) ) 

73 IF (ZZ.GE.ZO(N-l) ) GO TO 75 
N - N -1 

GO TO 73 

75 R = (ZZ -ZO(^-i) )/(20(N) -ZO(N-l)) 

E » «+EC(N) +(1. -R)+fcO(N-l) 

77 M - NX +2 -X 

G(I,(<Y + 1,L) ■ G(M,KY-1,L) -t 
G(M,KY+1,L) - G(2,KY-1,L) +E 

GK2(M,KY) - GK1(M,KY) 

GKl(f',KY) » G(M,KY,L) 

G(N,KY,L) - G(I,KY,L) 

GO TO 71 
END 
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SU8RCUTlNt Xi,wtfP 
CGlUrN RELAX4TICN 

CGMiGN G (192»^6,'♦)^ S0a93,3L»)i E0(131)i Z0(131) ^ 

IV(193>3f>), n rl(3&) , ITc2( 35 )# 
A0(193),Ai(i93)^A2(193),A3(193)j 
8C(25)#Bi(2o),B2(26)/B3(26)^ 
Z(35)#C1(35)#C2(35) ,C3(35)> 

XC(35),XZ(35),XZZ(35)#rC(35)»YZ(35),YZZ(35)j> 
NX# NY^NZ^KTLlf KTf 2# 1 S YM# K S Y^l# SC AL# SC AL Z # 
YAw#CYAW# SYAw# ALPHA#CA# SA#FMACH#N1#N2# N3# IG 
COMMGN/FLO/ S 1 P 1 » # Pi# P 2# P 3# B F TA # FR# I R# JR# KR# 0 6# I G# J G# K6# NS 
CGMMGN/SrtP/. GK1(193#26) #GK2U93#26)# 

SX(193)#S2a93>#SXX(193)#SXZ(193)#5ZZa93)# 

R0a93)#Ria93)»C(193)#0(193)# 

GIOJ 26) #G20(26)#G3O(2b) » GAC ( 26 ) » G1 ( 26) # G2 ( 2b ) # 
II# 12#K#L#N-)#LA#NX#KY»MY# Tl» AAO# 01# C2# T YAW# SI 



N 

a 

siO 


JI 

3 

c 


IF (FNACH.GE 

.1.) J1 


C ( Jl- 

1) 

0, 


D( Ji- 

i) 

G. 


S 

s 

1. 


11 

« 

1 


1 

DO 12 

s 

J»2#KY 

12 ti 


RC( J ) 

= 

1. 


Rl( J) 

■ 

1. 


Gl( J) 

S 

G10( J ) 

12 

G2( J) 

m 

G2C( J) 

21 

IP 

r 

1 411 


IM 

B 

I -ii 


J2 

a 

KY 


IF (IV( I#K).LT.2*AND.I.GT.LX) J2 » NY 
LV « lAFSd -1ABS( IV( I, H ) ) ) 

RG(KY) = AMlNodV# IABS( IV( i#K) ) ) 

Rl(KY) - LV 

AC « Tl + AKD+CKK) 

DU 32 J«J1#J2 

Ab = “Tl+Al ( I )*ai( J) 

BC « -Tl + BK J)^Cl(K) 

YP « SO(I#K) +aO(J) 

A « 1. -R0(0) +AG(I)»AO(I) 4YP+YP 

H » P 0 ( J ) / A 

FH » R0( J ) *A 

P = A0(1)4(A.*YP+YP -FH) 

C - YP4 (A.+AG ( l)+AO( I ) -FH) 

A = XZ(K)+XZ(!U -YZ(K)+YZ(K) 

B - (XZ(K) +XZ(K) )4YZ(K) 

AZ = -aO(I)+XZ(K) -YP+YZ(K) 

BZ » -AC(1)4Y2(R) +YP+XZ(K) 

CZ = H4H4(P+A -Q4B) -AC(1)4XZZ(K) -YP»YZZ(K) 

DZ « H*H4(Q+A +P+B) -AO(I)*YZZ(K) +YP+XZZ(K) 

DGI * S+(6{IP#J#L) -Gl(J)) 

DGJ » G(I#J + 1»L) -G(I#J-1#L) 
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r -4 tVj CO 


DGK 

DGII 

1 

DGJJ 

1 

DGKK 

1 

UGiJ 

i 

DGIK 

1 

DGJK 

1 

GX 

GY 

U 

V 

w 

1 

AU 
AV 
CXY 
QO 
AA 
HZ 
FXX 
FYY 
- FXY 
BV 
UU 
VV 

ww 

uv 


G(I>J#L+1) -GK1(I#J) 

G(I+1^J»L) ”G(i»J»L) 


+A3(i)*DGl 

G(l#J + lfL) —GCl^J^L) ”*G(I#J»L) 
-B3( J )*DGJ 

G(I>J,L + 1) -G(i,J,U -G(1>J#L) 
+C3(K )^DGK 

G(I + i>»j + l#L) “G ( i“ J ♦ 1 ^ L ) 
-G(I+1, J-i, L) +G(I-i, J-i,L) 
G(I+1.J,L+1) -6(1+1, 


♦ 6(1-1, ,J,L) 
+ G(i,4-l,U 
+G(1, J,L-1) 


-G(I-1, J,L+1) +6(1-1, j,L-l ) 

6(I,J+i,L + i) -G(l,J-i,L + U 
-G(l, J+1,L-1) +GU,J-l,l-n 
Aid )40G1 
-BK j)*DGJ 

GX -SX( 1)+GY +CA+AO( I) +SA+YP 

6Y +SA*AO(l) -CA*YP 

PO( J)^(C1(K)+DGK -bZ(I)+GY +SYAW 

♦CA+XZ(K) +SA+YZ(X) +H+(U»AZ +V+BZ)) 


U +,* + AZ 
^ +W^BZ 


H+(L+U +V*V) 

QXY +W+W 

oim(aag, . z+at) 

AZ+SX(1) -BZ +FH+SZ(I) 

1 . ♦H^AZ’^'AZ 

1. +SX(I)+SX(I) +H+HZ+hZ 
SX(1) +H+AZ+HZ 
AV -AU*SX(I) -FH+W^SZd) 


H^AU+AU 


H+BV+BV 


FHAW+W 


H+AU+BV 


UW ■ AU+W 

vw - av+w 

AXX ■ Rl( J ) + (FXX>AAA -fU) 

AZZ « FM+AA -WW 

AXZ - (R 1 (J) +Ri( J) )+( AZ*AA -UW} 

R » -(AXX + SXX( 1 ) +AZZ+SZZ(I) + A X Z* S XZ ( 1 ) ) 4 G Y 

+ Tl+( AA + (CZ+GX +(DZ -SX (1 )+CZ)+GY) 

-H+(CA+( AU+Ab -AV+AV) +(SA +SA)+AU+AV 
-aXY*(U+AC(I) +V+YP 
+(W +W)+( AO(I)+AZ +YP»BZ))) 

-WW+ (CA^XZZ (X) +SA + YZZ(K)) -W + ,i'»(U+CZ +V + DZ)) 
AXT • AeS(AU*Al(I ) ) 

AYT * ABS( 6 V*B 1 ( J) ) 

AZT ■ ABS(Fri+W*Cl(K) ) 

A = PO( J)*BB 1 A+AA/AMAX 1 ( AX 1 ,AYT,AZT, ( 1 . -RC(J))) 

AXT « A+AXT 

AYT « A+AYT 

AZT « A+AZT 

IF (OQ.GE.AA) GO TO 33 
AXX « AXX-»A 2 (I) 

AYY « (FYY+AA -VV)+B 2 (J) 
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I 


1<r-i 


t^ll 

AXY 

AXZ 

AYZ 

BP 

BK 

B 

R 

1 

6D R 3 5 
33 NS 
AXX 
AYY 
AZZ 
AXY 
AXZ 
AYZ 
bXX 
BYY 
BZZ 
BXY 
BXZ 
0YZ 

ALi 

PELTAG 

1 

D&II 

1 

DGJ J 

1 

DGKX 

1 

DGIJ 

1 

DGIK 

1 

DGJK 


- 4ZZ*C2(K) 

= -KK J)^(FXY*AA +LV)+(AB +A3) 

» AXZ*AC 

« -PK J) + (hZ*AA +VW) + (BC +BC) 

■ AYY 
= AYY 

» -AYY -AYY -Q1^(AXX +AZZ) 

« AXX*OGiI +AYY+DGJJ +AZZ+DGKK 

♦AXY+OblJ ♦AYZACGJK +AXZ*0G1K +R 

» NS +l 

■ UU+A2(1) 

« VV+bZ(J) 

* WW*C2(K) 

= b.’fS + Uv + AB 

« 8.+S+UW+AC 

■ B.+VW+BC 

- (FXX+OU -jU)*A2(1) 

= (FYY+00 -VV)»32{J) 

« <FH*CQ -»iW)AC2(K) 

«= -(FXY+QG +UV)+(AB +AB) 

- <AZ+00 -Uw)*(aC +aC) 

= -(HZ + OQ -*-Y*l) + (bC +BC) 

« AA/tQ 

» BXX + DGII 4-BYY + DGjJ +BZZ»CGKK 
♦BXY+OGl J 46 YZ»DGjK +3XZ406IK 

* G(UJ»L) -G(IM,JR) -G(R»J,L» +G2(U) 

+A3(I)+DG1 

■ G(1#J#L) -G(l>J-i#L) +G(I>J-2>L) 

-B3( J )*DGj 

* G(l,o,L) -G(I,J,L-1) -G{I,J#L-1) +GK2(I#J) 

+ C3{K)»DGi<l 

= G(I, J»L) -G(IM, J^L) 

-Gd^J-l^L) 

- G(I,JR) -6(I,J#L-1) 

-GUf^jJfL) ■►G( IP, J,L-1) 

« G(I, J,U -GlI, J,L-I) 


1 

GSS 

1 

BP 

Bh 

B 

1 


-G(1,J-1,L) +G( 1, J-l, L-i ) 

AXX^DGII +AYY-^DGJJ +AZZ>t'DGKK 
+AXY+D6IJ +AYZ+LGJK +AXZ*OGiK 
AQ+BYY 

BP -(AQ -l.)*(AYY +AYY +AXY 
-BP -BP -Q2^AC4(EXX +BZZ) 

+(AC -1. )+(2.^( aXX +AYY +AZZ) 


K 

35 IF 
FR 
IR 
JR 
KR 

37 R 
1 

B 

BM 


(A0S(R) 


(At -i.)»GSS 
LE.AES(FR)) GO 
R 
I 
J 
K 
K 


+ AQ>*L'tLTAG 
TD 37 


+ R 


B 


-AXT^CGKJ) -G(iM,J,D) 
-AZT’^IGKK I, J) -G{l,J,L-i)) 
-AXT -AYl -AZY 
+ AYT 


+ AYZ) 
+ AXY 


+AYZ +AXZ) 
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“B«*C ( J-1) ) 


(A8S(CG) 


B 

C( J) 

32 D(J) 

C6 

J 

DG M M-Jl 
CG 
IF 
DG 
16 
JG 
KG 

G2(J) 

GKJ) 

GK2(1# J ) 
GM( I, J ) 
G(Xf Jf L) 

^2 J 

IF (1V(I,K) 
A 
H 

FH 
AZ 
BZ 
HZ 
Fvy 

FXY 

DGI' 

DGK 
V 
U 
W 

G( I»KY+1»L) 


l./O 
e*ep 
B + (F 

0 . 

J2 

J2 

D( J ) 


-BM+0( J-1) ) 


-C( J)+CG 

GO 


TO ^3 


-CG 


T3 


51 

>A0( i)*A0( 1 ) 


4-SC ( I»K) *S0 ( I»K) 


•LF.ABS(DG) ) 

CG 
1 
J 
K 

Gl( J) 

G(I> J>U 
GKK J) 

G(I> J>L) 

G( I# J^L) 

J -1 
.LT.2) 60 
i, -kO(KY) 

R0(KY)/A 

KC(KY)»A 

-A0(1)+XZ(K) -SC(I/K)*YZ{K) 

-A0(I)+YZ(K) +S0(I»K)«XZ (K ) 

AZ+SXCl) -BZ +FH4SZ(1) 
i. 4SX(i)+5X(i) +H*HZ*HZ 
$X(1) +H+AZ4HZ 
S+(G{ IPiKYf L) -G2(KY)) 

G(I»KY,L+1) -GK2(1#KY) 

SA+AC(I) -CA4S0(I/K) 

A1(1)*DGI +CA4AG(I) +SA4S0(1^K) 

C1(K)4DGK +SYAW +CA+XZ(K) +SA4YZ(K) 

«= G(I,KY-1#L) 

•MV+d. -H + BZ4HZ) -UAFXY -W + HZ ) /( FYY^BU KY ) ) 


IF ( I.NE.ITEKK) ) GO TO 61 


M » NX +2 -I 

t ■ G(M>KY>L) -G(i>KY#L) 

NO » NO +1 

EO(NC) » EO(NC) 4P34(E -bO(NO) 

N « NO 

GO TO 61 

51 IF (l.GT.LX) GO TO 61 
E » 0 . 

IF (I V( I>K).N£.l) GO TO 57 

11 » Z(K) -TYa^4(XC(K) +S14A0( I)»A0( 1) ) 

53 IF (ZZ.GE.ZO(N-l) ) GO TO 55 
N ■ N -1 

GO TO 53 

55 R ■ (ZZ -ZO(N-l) ) /(ZO(N) -ZO(N-D) 

E * R + EC(N) -Mi. -R)4F.Q(N-i) 

57 M - NX +2 -1 

G(I»KY41,L) - G(M>KY-1»L) -B 
G(MiKY>lM) ■ G(l>KY-lfL) +1 
6K2(f^>KY) « 6K1(N,KY) 
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GKKN^KY) = 6(M#KY#L) 

GCM^HYiL) - G(IiKY^L) 

ei IF (I.EQ.NX) GC TG 71 
IF (I.6Q.2) RETURN 

I » I +II 

GU TC 21 

71 S « -1. 

II = “1 

I ■ II -1 

00 72 J*2#KY 

Gi( J) ■ G30( J) 

72 G2( J) * G^C( J) 

GO TO 21 

END 


+ E 


SUBROUTINE 

CALCULATES 

COMMON 


1 

2 

3 


7 


VELG (K»L^SV^SM>CP»X>Y) 

SURFACE VELOCITY 

G(193»2fc/4)»S0(19 3,35),E0a31)#Z0(13i)> 

IV(193,3S)»ITtU3f;)^lTF2(3E)# 

AC(193)»a1(1‘»3)>A2(193),A3(193), 

ao(26)fei(26)»B2(26),B3(26)> 

Z(35),CH35),C2(3S)/C3(3iJ)# 

XC(35),XZ(3^)#XZZ(35),YC{35)#YZ{35),YZZ(35)j 
NX, NYf NZ»KT£1» KTE2» iSYh, KSYM# SCAL> SCALZ# 

yaw,cyaw,syaw»alfha,ca, sa,fmach,ni,N2,n3, id 


DIMENSION 


SV(1),S^(1) 

,cpa),x(i),Y(i) 

11 

a 

ITEKK) 


12 

a 

ITE2(K) 


J 

a 

NY 41 


01 

= 

.2+FMACH442 


T1 

a 

l./( .7*FMACH*+2) 

DU 12 1*11, 

12 


FH 

H 

IF (IV(I, 

a 

AC(1 )4A0( I ) 

0 • 

. N E . 0 ) H » i 

+S0( l,K)+50( I,K) 

K) 

./FH 

AZ 

a 

-AC(I)4XZ(K) ~SOU,K) + YZ(K ) 

BZ 

a 

-ACm+Y2(K) 4S0U,K) 4aZ (K ) 

DSI 

a 

SC(I+1,K) 

-S0( 

DSK 

a 

S0{I,K+1> 

-S0U,K-1) 

SX 

a 

A1 { 1 )»DSi 


sz 

a 

C1{N)+USK 


DGI 

a 

GU + I, J,L) 

“•G ( 1 — i, »J, L ) 

DGJ 

a 

G (1, J+1,L) 

“GU, J“1,U 

DGK 

a 

G( I, J,L+1) 

-G(I, J,L-1> 

U 

a 

A1(I)40GI 

♦ SX481 ( J)*DGJ 4CA4A0(I) -fSA^SOn 

V 

a 

-ei( J)40GJ 

+SA*At(I) -CA4S0(X,K) 

W 

a 

CKK )+DGK 

+ SZ4BK J)+OGJ +SYAW 



+CA+XZ(K) 

+SA4YZ(K) +H*(U+AZ +V*BZ) 

00 

a 

H+(U’*U 4V + \/) tW^'W 

SV( I) 

a 

SIGN{SQRT(QQ),U) 

IF { IV( I, 

K) 

. bo.o) sv( n 

- SV(I-l) +SVU-1) -SVd-2) 
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QQ 

SM( I) 
CP(I) 

xm 
12 Y(I) 
RETURN 
END 


1. +Q1*(1. -QQ) 

FMACH+SV< I)/SQKT(QQ) 

Tl+(QQ + t3.i? -1.) 

XC(K) +,3 + SCAL»(A0<I)'*‘A0(I) -S0(I#K) ♦SOC I»K) ) 

YC(K) ♦SLAL*A0n) + S0(l»K) 


SUBROUUNb CPLCT ( 1 1 » U# F MAC H, X , Y, C P ) 

C PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 

DIMENSION KQ0E(2) »LlNE(lcC),X(l),Y(l)»CP(i) 

DATA KODt/lH ^1H+/ 

IwRIT - 6 

WRITE (IwRlT»2) 

2 FDRMAK 5OE0PL0T OF CP AT cQUAL INTcRVAlS In THE MAPPED PlANE/ 
1 lOHO X >10H Y ,1CH CP ) 

CPO * ((1. +.2*FMACH»*2) ♦*3.5 -1. ) / ( .7'^FMACH*+2) 

DO 12 I«l#100 
12 LINE(I) * KODE(l) 

DO 22 I-I1#I2 

K » 30. ♦(CPO -CP( I) ) 

K • MINCdOO^K) 

LINE(K) « KGD£(2) 

WRITF (IWR1T»61C) X(I)»Y(I)KP(l)#UNt 
22 LrNE(K) » KQDE(l) 

RETURN 

610 FORMAK 3Fi0.4#lC0Al) 

END 


SUBROUTINE FORCE ( II, 1 2^ X, Y, CP , AL^ ChQR 0» XM, CL> C 0# CM ) 
C CALCULATES SECTION FORCE CQcFFICIEnTS 

dimension X(1), Y(1),CP(1) 

RAD * 57.295779^130823 

ALPHA » AL/RAO 

CL = C. 

CD - 0. 

CM “0. 

N » 12 -1 

DO 12 I=Ii,N 

DX ■ (X(l+1) -X(1))/CHQRD 

DY » (Y(I+1) -Y(1))/CHQRD 

XA - (.5*{X(1+1) ♦X(I)) -XM)/CnaRD 

YA = .5+(Y(I+l) +Y(I))/CHOKD 

CPA - .5 + (CP(I + l) +CP(I) ) 

DCL * -CPA+DX 

DCO » CPA+DY 

CL » CL +DCL 

CD ■ ro +DCD 
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CM 

* CM +DCD+YA 

-OCL + XA 

DCL 

= CL^CUS(ALPHA) 

-CD'< SiN( ALPHA) 

CD 

- CL^SIN( ALPHA) 

♦CD«CDS( ALPHA) 

CL 

RETURN 

tNO 

- OCL 



c 


SUBROUTINE TOTFOR (KTElf KTt£> CHUhC»SCL> SCO, SCM,Z^XCf 
1 CLfCD^CMP^Cf'R^CMY) 

CALCULATES TOTAL FORCE CCEFFlCltNTS 

OnENSION CHORD( i), SCL(l)f SCD(1 ), SCM( i),Z( i),XC( 1) 
SPAN = Z(KTE2) -Z(KTEi) 

CL - C . 

CD '0. 

CMP » 0. 

CMR = L. 


ChY ■ 0. 

S » 0. 

N » KT£2 -1 


DO 12 K=KTE1,N 

DZ » . ( Z{K + 1) -Z (K) ) 

AZ =* .S>i^tZ(K + l) +Z(N)) 

CL » CL +DZ>» { sen K+l) YCHuRD( ) +SCL ( K )»CHQRD (K) ) 

CD = CO +DZ^( SCD(K + i) ♦CHGRD( K + i ) +SCD ( K ) *CHORD (K. ) ) 

CMP ■ CMP +DZ+(CHCKD(K + l)»(SCM(K + l)*CrlQRD(K + l) 

-SCL{K+i)+XC{K+l) ) 

♦CHORD ( K) *( SCM(K)+CHORO(K) 

-SCL(K)+XC(K) ) ) 

CMR = CMk +aZ4DZ»(SCL(K+1)+CHURD(K+1) +SCL ( K) ♦CHORD (K ) ) 

CMY » CMY +A7 + uZ+(SCD(K + l)+CHQRC(K + l) + SCO ( K ) ♦CHORD ( K ) ) 

S = S +DZ»{CH0RU(K+1) +CHGkD(K>) 

CL - CL/S 

CD = CD/S 

CMP - CMP*SPAN/i**2 

CMR » (CMR HCMk) / ( S + SP AN ) 

CMY » (CMY ♦CMY) /(S^SPAN) 

return 

END 


SUBRCUTINE 
HALVES MESH 
COMMON 


R EFiN 
SIZE 

G{193>26fA),SO(i9B^3S)^eO(131)#Za(131)^ 

lV(193,3S),ilEl(35),ITfc2(3S)f 

A0(193)^Al(193),A2(l93)»A3(193j> 

60(26)^81(2 6 )»82(26)>B3(26)t' 
Z(35)>C1(3^),C2(35),C3(35), y 

XC05), XZ(35)>XZZ (3S)» YC ( 3S) »YZ( 35), YZZ(^)» 
NX,NY,NZ,KT61,KTE2, 1SYM,KSYM»SCAL,SCAL^ 
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7 YAW^ CYAW^ SYAW, ALPH^ ,CA, SA,FMACH, Nl#N2f N3> ID 


MX 

S 

NX +1 


KY 

■ 

NY +1 


MY 

B 

NY +2 


MZ 

m 

NZ +3 


MXO 

- 

NXZ2 

+ 1 

MYO 

m 

NY/2 

+ 2 

MZO 

m 

NZ/2 

♦ 1 

K 

m 

i 


IF (KSYM. 

EQ 

.0) GO 

TO 

MZO 

S 

NZ/2 

+ 3 


BUFFER IN (NX,1) ( G ( 1, 1 ) # G ( MX 0# M YO# 1 ) ) 
IF (UMT(Nl) .GT.O.) Gil TO 4C1 
K - 2 

11 BUFFER IN (Nl,l) (G(l»l>l)»G(hXCiMYO,l) ) 
IF (UNIT(Nl) .GT,0. ) GO TO ACi 
J «= NY/Z +1 

Jj » KY 

21 i « MXO 

II » MX 

31 G(II>JJ#i) = G(i»J,l) 

I -1-1 

II - II -2 

IF (I.GT.D) GO 10 31 

J - J -1 

j j * J J — 2 

IF (J.GT.O) GG TO 21 

DO'42 J»1#KY>2 
00 AZ I*2#NX#2 

42 G(I^J»1) - .5*(G<I+1> J>i) +G(l~l,J,i)) 

DO 52 I-1#MX 
DO 54 J»2>NY>2 

54 G(I#J»1) - .5+(G(i# J + l#l) +G(I#J-l,D) 

52 G(I,hY»l> » 0, 

BUFFER OUKNZ.l) ( G ( i, x, i ) , G ( MX # MY, i ) ) 

IF (UNIT(N2) .GT.O. ) GO TO 4C1 
K - K +1 

IF (K.LE.MZO) 60 TO 11 
REWIND Ni 
REWIND N2 

BUFFER IN (N2,l) ( G ( 1 , 1, 1 ) , G ( MX , MY, 1 ) ) 

IF (UNIT(N2).GT.O.) GO TO 401 

BUFFER IN (N2,l) ( G( 1, 1, 3 ) , G( MX ,MY, 3 ) ) 

IF (UNIT(N2) .GT.O. ) GO TO 401 

BUFFER 0UT(N1,1) ( G( 1 » 1, 1 ) , G ( MX , MY, 1 ) ) 

IF (DNIT(Nl) .GT.O. ) GO TO 401 
K -1 

IF (KSYM.NE.O) K - 2 

111 K » K -*•! 

DO 112 J»1,MY 

DO 112 1=1, MX 

112 G(I,J,2) - .5*(G{I,J,1) +G(I,J,3)) 

DO 122 L=2,3 

BUFFER OUT(Ni,l) ( 6 ( 1, I, L ) , G ( MX, MY, L ) ) 
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IP (L'NIT(M).GT.O.) GO TO 401 
122 CGNTiNUE 

IF (K Gf rc 2Ui 
DU 132 J«1,MY 
DO 132 i=l^MX 
132 GdiJ/l) « G(1#J#3) 

EUFFtR IN (H2»l) ( G ( 1 # 3 ) # G ( MX # MYi 3 ) ) 
IF (LMT( N2) .GT.O. ) GO TO 40i 
GO TC Hi 
201 RfcWINO M 
REWIND N2 
DO 202 L-l#3 

bUFFER IN (Nl»l) (G(l»i>L)#G(hX,MYH) ) 
IF ( UMT(M) .GT.O. ) Gj TO ‘tOl 
2C2 CONTINUE 

BUFFER 0UT(N2H) ( G ( IHH ) # G ( MX , MY» i ) ) 


IF (LNIT(N2>.GT.O.) GO TO 401 


T Y AW 

ii 

NO 

EC (NO 
K 


(K SYM.Nt 


IF 

211 N 

I 

IP 

II 
12 

DO 
DSI 
DSK 

sx 

bZ 

R 

A 

H 

FH 

AZ » 

BZ 

HZ 

FYY 

FXY 

DGI 

DGK 

V 

U 

W • 

212 G(I>KY+1,2) 


CO TQ 251 


GO TO 231 


-SD(l-i#K) 

-SO(i>K-l) 


bYAV ^CYAw 
• i»»SCAL 
KTEl -1 

0 . 

2 

0 ) 

■ NO 

» MXO +1 

(K.LT.ME1.0P.K.GT.KTE2) 

■ ITEI(K) 

- ITE2(K) 

212 l-ll»I2 

= S0(I+1>K) 

■ S0(I»K+1) 

» Al(l)+Dbi 
» C1(K)+DSK 
= A«lNC(iHV(I,K) ) 

• 1. -R +AO(i)+AG{J) 

« R/A 

■ R*A 

» -A0(1)+XZ(K) “b( (1^K)+YZ(K) 

• -AO(I)+YZ(K) +SO(HK)*XZ(K) 

« AZ+SX -ez +Fh*SZ 

1. -PbX + bX ♦H*HZ*HZ 
SX -RH^AZ + HZ 

G(I+ifKY,2) -G(I“1»KY,2) 
G(I,KY»3) -G(I»KY>1) 

SA+AC(I) ~CA+SG(1#K) 

Aim + DGl +CA*A0(1) -RiA^SUlHK) 


+SC( I^K)»SO( I#K) 


NO 

EC(Nt ) 

N 

I 


C1(K)+DGK +SYAW +CA+XZ(K) +SA*YZ(K) 

« G(HKY-1»2) 

4(V+(1. -H+BZ4HZ) -U+FXY -W^HZ ) /( FYY+Bl ( NY) ) 
NO +1 

G(12»KY>2) -G(I1>KY,2) 

NO 

II 


IF (K»NE«KTE2»DF»YAV/»Lc»0») GO TQ 231 
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221 I - I +1 

M « NX +2 -I 

NO « NO +1 

EO(NC) » G(«#KY»2) -0n*KY,2) 

IF (I.LT.MXO) 60 TO 221 
i » II 

231 I » I -1 

£ - 0 . 

IF ( 1V( r,K) .NE.i) GO TO 237 

ZZ = Z(K) -TUw*(XC(K) +S1^A0( I) ♦A0( 1 ) ) 

233 IF (2Z,GE.Z0(N-1) ) GO TO 23t? 

N « N -1 

GO TO 233 

235 ft « (ZZ -ZO(N-l) ) /(20(N) -ZO(N-i)) 

£ * RttC(N) +(1. -P)^£0(N-1) 

237 M • NX +2 -1 

G(I,KY+1#2) * G(M#KY-1,2) -E 
G(M#KY+1#2) « G(J.>KY-1,2) ♦£ 

IF UVdf K) .NE.-l) GO TO 241 

G(I»KYf2) = .^^G(I,KY»1) + . 2 5+ ( G ( I , K Y» 3 ) +G(M#KY/3)) 

IF ( IV( I>K + 1) .LT.l) 

IGd^KYfZ) - .£ + Gd»KY»3) + . 2 5 ♦ ( G ( I , K Yd ) +G(M,KYd)) 
G('1dY/2) = GddYd) 

Gd>KY-l,2) ■ .i*lC(I,KY,2) +G(dKY-2d)) 

G(?1iKY-ld) « .i^{G{M>KYi2) +G(M#KY-2d)) 

241 IF d.GT.2) 60 TO 231 

251 K = K +1 

IF (K.EU.MZ) GO TO 261 
DO 252 J«1,MY 
DO 252 

Gd.jd) » G(Idd) 

252 Gd,J#2) = Gddd) 

BUFFER 0UT(N2d) ( G dd» 1 ) > G ( MX , rt Yd ) ) 

IF (UNIT(N2).GT.G.) GO TO 401 

BUFFER IN (Nl#l) (G(ldd)d(MX,MYd) ) 

IF (LMT(M) .GT.O.) GO TO 401 
GO TO 211 

261 EG(NC+1) = 0. 

DO 262 L»2d 

BUFFER 0UT(N2d) ( G ( 1 d# L ) > G ( MX # M Yd ) ) 

IF (UNIT(N2) .GT.O. ) GO TO 401 

262 CONTINUE 
REWIND N1 
REWIND N2 

00 3C2 K»1#MZ 

BUFFER IN (N2d) ( G ( 1# Id ) d ( MX, MY# 1 ) ) 

IF (UNIT(N2) .GT.O.) GO TO 401 

BUFFER OUT(Nid) ( G ( 1 d, 1 ) , 6 ( MX , MYd ) ) 

IF (UNIT(Nl) cGT.O.) GO TO 401 
302 CONTINUE 

IC = 1 

RETURN 

401 10 = 0 

RETURN 
END 
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c 


SU^RCLTCNfc 5MQ0 


1 

2 

3 

u 

b 


SMOOTHS POTtNTlAL 

COMMON G{19 3^2 6»'t)#S0(193»35),E0(131)>ZU(131)i 

IV(193»35),ITEil3i>),ITE2(35), 
Au(193)>Ai(i93),A2(193)#A3(193), 
B0(25),tl(26),b2(26),B3(26), 
Z(35)^C1(35),C2(35)^C3(35), 
XC(3^),X2(35),XZZ(35)#YC(35)#YZ(35)^YZZ(35) 


t 


MX 

S 

NXpNYpNZpKTEIpKTcZp ISYMpKSYMp SCALpSCALZp 
YAW, C YAW, S YAWP ALPHA pCAp SApFMACHp N ip N2p N3p 10 
NX +1 

KY 

= 

NY 

+ x 

MY 

M 

NY 

+ 2 

MZ 

m 

NZ 

+ 3 

K1 

M 

2 


K2 

a 

NZ 


IF 

(MSYM.fcO 

.G ) 

GO TO 1 

K1 

u 

3 


K2 

c 

NZ 

+ 2 


1 PX « L./6-. 

PY » i./t. 

PZ 

DO 2 L-l»3 

BUFFER IN CUfl) ( 0 ( i # L » L ) > G ( MX # MYi L ) ) 
IF (LMT(M) .GT.O. ) 60 TO 51 

2 CONTINUt 

eUFFbP aUT(N2#l) <G(1»1#1),G(MX,MY^ 1) ) 
IF (LNIT(N2).GT.U.) GO TO 51 
K ■ K1 

11 K - K +1 


DO 12 J=3,NY 


DO lA 1-2, NX 


G(I,J,A) 

1 

2 

3 

G(l, 


(1. -PX -PY -PZ)+G(1»J,2) 

♦ •5^PX+(G{ 1 + 1, J,2 ) +G(I-1,J,2)) 

+ .5+PY+(G{ I, J + i,2 ) +G(I,J-i,2)) 

+ .5^PZ^(G(i, J,3) +G(I>J,D) 

G(l# J#2) 


12 G(MX#J^^) - G(MX,J#2) 
DO U I=1>«X 


G(I#1#A) - G(U1#2) 


G( Ip2» A ) = G( Ip 2»2) 

G(IpKYpA) - G(IpKYp2) 

16 G(I,MY,4) = G(I,MY, 2) 

BUFFtK aUT(N2,l) ( G ( 1, 1, 4 ) , G ( MX , MY, 4 ) ) 
IF U NIT(N2) .GT.O.) GO TO 51 
IF (K.FU.K2) GO TO 31 
DO 22 J-IpMY 
DC 22 1=1pMX 


G(IpJp.I) ■ G(1,J,2) 


22 GII,J,2) = G(IpJp3) 

eUFFFR IN (Nlpl) (G(1 p1p 3),G(MX,MY,3) ) 
IF (LNIKM) .GT.O. ) GO TO 51 
GO TC ii 

31 BUFFER QUT(N2p1) ( G ( 1 , 1, 3 ) p G ( MX , f1 Yp 3 ) ) 
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o o r> 



IF (UNITCN2) .GT.O.) GC TO 51 
REWIND HI 
REWIND N2 
DO 42 K-1#MZ 

BUFFER IN (N2#l) ( G ( i , i ) , G ( MX i ) ) 

IF (UMT(N2) .GT.O.) GO TO 51 
BUFFER OUKNi^l) ( G ( 1 , i, 1 ) , G ( MX , MY, L ) ) 
IF {LNIT(NI) .GT.O. ) GO TO 51- 
42 CONTINUE 

10 “1 

RETURN 

51 10 » 0 

RETURN 
END 


SUBROUTINE SPLIF(M,N,S,F»FR,FRP,FPPP,KM,VM,KN, VN, MCUE,FUM>IND) 
SPLINE FIT - JAMESON 

INTEORM PLACED IN FPPP IF MODt. GREATER THAN C 

INO SET TO ZERO IF DATA ILLEGAL 

DIMENSION S(l),Fn),FP(i),FPP(l)#FPPP(i) 



IND 

« C 


K 

« iA6S(N -M) 


IF 

(N -1) 81,61,1 

1 

K 

= (N -M)/K 


I 

« M 


J 

» iM +K' 


DS 

■ S(J) -SU) 


D 

* OS 


IF 

(DS) 11,81,11 

11 

DF 

• (F(J) -F(I))/DS 


IF 

(KM -2) 12,13,14 

12 

U 

* . 5 


V 

» 3. ♦(OF -VM)/OS 


GO 

TO 25 

13 

U 

■ 0 . 


V 

» VM 


GO 

TO 25 

14 

U 

■ “1 . 


V 

= -DS^UM 


GO 

TO 25 

21 

I 

■ j 


J 

= J +K 


DS 

» S ( J ) -S( I ) 


IF 

(D^OS) 81,81,23 

23 

DF 

- ( F( J) -F(l) )/DS 


B 

« l./(DS +DS ♦U) 


U 

* B +DS 


V 

- 3+(fc.+DF -V) 


25 FP(I) * U 

FPPd ) -V 

U » (2. -U)*DS 
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V 

= 

6.+LF +DS^V 


IF ( J - 

N) 

21,31,21 

31 

IF (KN 

-2) 

32,53,3A 

32 

V 

GO TC 35 

3 

(6.+VN -V)/U 

33 

V 

GO TO 35 

C 

VN 

3A 

V 

S 

(DS+VN +FPP(I))/(1, +FP(I)) 

35 

B 

m 

V 


D 

m 

DS 

A1 

DS 

s 

S(J) -S(I) 


U 

m 

FPP(l) -FP(X)+V 


FPPP( I ) 

c 

(V -U)/ES 


FPP( I ) 

m 

U 


FP(I) 

= 

(F(J) -F(I))/DS -DS*(^ +U 


V 

■ 

U 


J 

s 

I 


I 

m 

I -K 


IF U - 

M) 

^l,51,Ai. 

51 

I 

X 

N -K 


F H P P ( N ) 

a 

FPPP( I ) 


FPP(N ) 

a 

3 


FP(N) 

a 

DF ♦D*(FPP(1) +C ♦8)/6. 


INO 

a 

1 


IF (MODE 

) 81,31/61 

61 

FFPPl J ) 

a 

FQM 


V 

a 

FPPl J) 

71 

J 

a 

J 


J 

a 

J +K 


DS 

a 

3(J) -S(l) 


U 

a 

FPPl J ) 


FPP? ( J ) 

a 

FPPP(I) +.5+0S^(F(I) tFlu) 


V 

a 

U 

81 

IF < J - 

RETURN 

END 

N) 

71, 61,71 


+U) /6. 


DS+Di*( U 


+V> /12.) 


SL'SRLUTiNfc 1NTPL(H1,NI,SI>F1,M,N, F,FP,FPP,FPPP, MODE) 

C INTEFPQLATIQN USING lAfLDR SERIES - JAMESON 

C ADDS CORRECTION FDR PIECtWISF CONSTANT FOURTH DERIVIATIVt 

C IF MODE GREATER THAN 0 

DIMEhSION Sia),Fl(l),SU)#F( 1),FP (1) ,FPP{1), FPPP(l) 

K = IABS(N ~M) 

K « (N -M)/N 

I “ M 

MIN » MI 

NIN * N1 

D - S(N) -S(M) 

IF (r+{SnNl) -SKMI))) 11,13^13 
II MIN - NI 

MN « Ml 
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13 

KI 

» IABS(NIN -MIN) 


IF (Kl) 

21,21,15 

15 

KI 

» (NIN -M1N)/Kl 

21 

11 

■ MIN -KI 


C 

» 0. 


IF (RODE) 31,31,23 

23 

C 

■ 1 • 

31 

11 

* II +KI 


SS 

■ SKID 

33 

I 

•I +K 


IF (I 

-N) 35,37,35 

35 

IF (D»(Sd) -SS)) 32, 33, 37 

37 

J 

= I 


1 

» I -K 


SS 

- SS -Sd) 


FPPPF 

- C + (FPPP(J) -FPPPd ) )/(S{ j) - 


FF 

- FPPPd) +. ZS + SS^FPPPP 


FF 

« FPPd) +SS*FF/3. 


FF 

« FPCl) +.:/tS5 + FF 


FKII) 

■ Fd) +SS + FF 


IF ( 11 

-MN) 3i,Al,3i 

A1 

RETURN 

END 



S( i ) ) 


SUBROUTINE RPLOT ( I PlUT> NP h S > RfcS^ C GUN f ^ T1 TL FM AC H, Y A/ A l » 
1 Ni>N2#N3) 

C PLOTS CGNVERGEKCE RATE 

UIMENSION RESd )#CDUNT( 1), T ITl t ( 20 ) # K ( 20 ) 

IF (NReS.LE.l) RETURN 
IF (IPLOT.EQ.O) GO TO 11 

CALL PL0TS8L( 100C>2RHANT0NY JAPESON 1096GAK) 

CALL PL0T(1.25d .#-3) 

11 IPLQT = d 

RATE » ( ABS(RES(NRES)/KES(1) ) ) 

1 ♦♦(l./(CGUNT(NRES) -CDUnT(I))) 

ENCarE( B0#12^ R) TITLE 

12 F0RMaT(20AA) 

CALL SYH80L(1.> .5# . lA,R,G.»aC ) 

ENCODE (50dA:» R) FMACri^YA^AL 

lA F0RMAT(5HhACH ^ F9 .3# A)<» 5HYAw > F9 . 3# A XoH A L P H A» F 9. 3 ) 

CALL SYMBDL(l.f.2:>^.lA#R»0.^5C) 

ENCQf e ( 32d6> R) RES ( 1 ) > RES ( NRES ) 

16 FOR-IAT { ;?HRES1 > L 9 . 3 » AX, i?HR E SE >E9.3) 

CALL SYMBGL(1.>0.,.1A»R,0.>32) 

ENCODE ( 50, 18, R ) COUNT ( 1 ), COUNT ( NR t S ), R A T E 
18 F0RMAT(5HR0RK1, F9.2,AX,5HW0RK2, F9.2, AX, SrIRATh ,F9.A) 

CALL SYMB0L{i.,-.25,.lA,R,G.,^0) 

ENCODE ( 2A, 20, R ) N1,N2,N3 
20 F0RMAT(6HGRI0 ,IA,3H X ,IA,jH X , lA ) 

CALL SYMBOLd .,-.5, .1A,R,0., 2A) 

RMIN = 0, 
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RHA< = 0. 

CUUNTl » CQUMd) 

kbSl « RbS(l) 

DO Zl I»lfNPeS 


CQUMTd ) 
R6S(I ) 
RMAK 
22 RMIN 
YSCAL 
YINT 


» CDUNT(I) -COUNTl 
• ALQG( ABS(R£S( I ) /kESl) } 

» AMAX1(R?1AX,R£S( I) ) 

- AMINI (RriN,Rf.S ( I) ) 

= i./AL06(i0.) 

■ 1 . 

IF ( YSCAL + K='1IN.LT.-6. ) Yif^T - 2, 

YLOW », -6.4Y1NT 

YSCAL « YSCAL/YIM 

XINF - 50. 

IF (CJUNT(NRFS) .Gr.3C0. ) XINT » IGO. 

IF (CCUNT(NPES) .GT.6C0. ) XINT » 200. 

IF (COUNT(NF£S) .GT.1200. ) XINT * 500. 
IF ( CQUNK NR6S) .GT.6000. ) XINT • 1000. 
XSCAl « I./XINT 


CAU PLQT( 

CALL AXIS(0.#-3.,10HLL3(£kkOF 1C> 8., 90.j YLOW# Y iNT# 0) 
CALL PLOT( 3.»-i. »-3) 

CALL AXiS(-3.,0.> AHNCYC>-4#b.#0.i0.» XINTjO) 

UC 32 I-1>NRES 


CDUNT(I) * XSCAL^CODNT(I) -3. 

32 RtS(J> • AN1M(2.,YSCALYRES( 4 ) ) +6. 

CALL LlNE{CQUNT»RtS?NKtS^i^0fl#0.#l.>0.»l. ) 
CALL PLOT ( 8. 5, -1.5, -3) 

RETURN 


LNO 


SUBROUTiNc GRAPH ( I F LuT, I 1, 1 2, X, Y, C P, T ITL E , F HAC H, Y A, AL , 
1 Z,CL, CD, CHDRDO, XSCAL, PSCAU 

C GENEFATES CALCOkP PLcTi 

DIMENSION X(l), Y(l) ,CP(1),T1TLE(20),R(20) 

IF (iPLOT.FQ.O) GC TO 11 

call PLOTSBU 10C0,2AHANTQNY JANESUN 1096C4R) 

CALL PLOT( 1.25,1 .,-3) 

11 IPLGT ■ 3 

ENC0DH80,12,R) TITLE 

12 F0RMAT(20A<^) 

Call symrgl( . 5,c.,.i^,k,o.,BL) 

ENCODE A4, is, R ) FI-ACH,YA,AL 

lA FLRMaTC 5HN ACH , F 7 . 3 , SX, 5 HY AW , F7 . 3, S X, 5H AL PHA, F 7. 3 ) 

CALL SYMBCK . 5, - . 25, . 1^, R, 0. , ) 

ENC0Dt(AS,16,P) Z,CL,CD 

It F0RMd(5HZ ,F7.2, AX, 5HCL , F7. A, AX, 5HC0 ,F7.4) 

CALL SYMPOK .5,-.5, .1S,R,0.,AA) 

XMAX • X(X1) 

XMIN - X(I1) 

YMIN « Y( II ) 
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.00 22 

XMAX « AhAXl(X(I),XMAX) 

XMIN - AM1M(X( i),XMIN) 

22 Yf«lN “ AMIM( Y< 1), Y«iN) 

SCALX ■ 5./<XMAX -X^^IN) 

IF (XSCAL.GT.O. ) SCALX » XSCAL/<XMAX -XMIK) 
IF (XSCAL.LT.O. ) SCALX «« AeS( XSCaU /CHG fiDO 


PINT 

= -.4 


IF (PSCAL.NE.O.) PINT • 

“A8S( PSCAL ) 

SCAlF 

» l./PINT 


PMIN 

» -3**PiNT 


PMAX 

» S.^J'PINT 


DO 2A 

I-Il, 12 


X( I) 

* scALx+(x(n 

-XMIN) +. 

24 Y(I) 

- SCALX^(Y(1) 

-YMN) ♦. 

CPMAX 

a 0» 


IMAX 

■ (12 +ID/2 


N 

« (12 -Il)/b 


M 

» XMAX -N 


N2 

» IMAX +N 


DO 26 

I«N1,N2 


IF (CP( I ) .LE.CPMAX) GQ 

TO 2^ 

CPMAX 

- CP(1) 


IMAX 

* I 



26 CONTINUE 

N -12 -ri ♦! 

CALL LlNE(X(Il),Y{Il)^iNfl,U,l,0.,l.,c-.#l.) 

CALI PLOKO., A.3>-3) 

Call AXIS( 0 .^- 3 .> 2 HCP# 2 > 8 ., 9 C.> PMN# PIM,C) 

CPC « (((5. +FMACH + + 2) /6. )* + 3. 5 -1 . ) / ( . T-^^F M AC H + + 2 ) 

IF (CPC.G6.PMAX) CALL S YN80L ( 0 . # SC AL P *C PC» . AO, 1 5#0 . #-l) 

DO 32 I«ri,INAX 

IF (CP( I) .LT.PNAX) go to 32 

CALL SYMbOL(X(I),SCALP+CP(l), .U7,3,A5.,“1) 

32 CONTINUE 

DO 3A I-IMAX, 12 

IF (CP(I) .LT.PMAX) 3U to 3A 

CALL SYMfiOL(X(I),SCALP^CP( i), .07,3,0.,-i) 

3A CONTINUE 

CALL PLOT(12.,-A.5,-3) 

RETURN 

END 


SUB'^LUTINE THREEDC IPLQT,SV,SM,CP, X, Y,TiTL£, YA, AL, 

1 VLO,CL,CO,CHCRDO,XSCAL,PSCAL) 

C GENERATES THREE DIMENSIONAL PLOTS 

COMMON G(193,26,A),S0(193,36),E0(131),Z0(i31), 

1 IV(193,35),irtl(3£),iTt2(35), 

2 A0(193),Al(i93),A2{i93),A3(193), 

3 B0(26), 81(26), B2<2E ),B3'26), 

4 2(35),Cl(35),C2(3!j),C3(35), 
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5 XC(3‘),XZ(35)#XZ2(35),YC(35)>YZ(35»>Y2Z(35)> 

6 NX,NY,NZ»KTti»KTt^>ISY^iKSYM.SCAL#SCALZ» 

7 YAW^CYA^# i YAk> ALPHA SA>fMACH»Nl#N2# N3f IQ 
DIMENSION X(1),Y(1),SVU)#SM( l),CP(l),ri7L£(2C)>R(20) 

LX - NX/2 ♦! 

MX - NX ♦! 

NY = NY +2 

IP ( XSC AL.NE. 0. > SCALX « . i + At».S ( XSCAL ) /CHQPOO 
if (PSCAL.GE.C. ) SCAlX « 5./(Z(KTE2) -Z(KT6D) 

SCALP '* -1.25 

IF < PSCAL.N6.C. ) SCALP = -.5/A8S ( PSCAL) 

SX -2. -SCALX+XC(KTE1) 

TX - 3.5 

IF (IPLQT.uC.C) GO TO i 

CALL ‘>LDTSBL( 10C0» 2AhAinTQNY JAMESON 109t0^R) 

CALL PLDT(1.25»I .,-3) 

1 IPLQT ■ 0 

N « 1 

ENCOD£{ 12»2#R) 

2 FQPMAT(12hVIEti OF niNG) 

CALL SYKdOL(2., .5».1A,R,0.#12) 

11 DO 12 L-l#3 

BUFFER IN (Nl,l) (G(l#i»L)»G(HX,MYf L) ) 

IF (CNIT( M) .GT.O. ) GO TO 101 

12 CONTINUE 

K =2 

21 K - K +1 

IF (K.GT.KTE2) C-Q TO fei 
DC 22 J»i>MY 
DO 2? I«1#MX 
GtI^J»l) ■ G(I.J>2) 

22 G( I, J,2) = G( I, J,3) 

eUFFI-R IN (Nl,l) (G(l,i,i),G(MX,MY,3) ) 

IF (iMT(Ni) .oT.O. ) GO TO 101 
IF (K.LT.KTEi) GO TO 21 

11 ■ I TE.1 (K ) 

12 = ITE2(K) 

CALL VELG (K>2f SV»SM,CP»A>Y) 

IF (K.6T.KTE1) GO TO ^1 
fcNC'UEI 3C#32>R) TITLE 
32 FDRr1AT(2CAA) 

CALL $YM8QL( .5,C .» .1A,R,0.#C0) 
bNCQCE(A4#3A»R) FflACrtfYA.AL 

34 FORMAT! 5HMACH > F 7 . 3 , 4 X , 5H Y Aw , F7. 3» A X, 5HALPHA, F7. 3 ) 

CALL SYMBOL! . 5»-.25,. 14,R#0.,44) 

ENCODE! 44#3e»R ) VLD#CL#CD 

36 FQRMAT(5HL/D • F 7 . 2 , 4X, 5HC L # F7.4> 4X> 5HCD >F7.4) 

CALL SYMeCL(.5f-.5» .14,R.0.f44) 

41 SY - 5.+ (2(K) -Z(KTED) /(Z(KTE2) -Z(KTED) +2.75 

DO 42 I-n,I2 

X(I) « SCALX+Xd) +SX 

Y( n » SCALX + Yd) +SY 

42 CPd) « SCALP + CP(d +SY 

IF (M.EQ.2) GO TO 51 
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N » 12 -II +1 

CALL LINc( X( II ), Y( ID ,0., I. ) 

GO TC 21 

51 N « 12 -LX +1 

CALL LINe(X(LX)#CP(LX),N,i#G»l»O.D.#O.D.) 
N - LX -li +1 

DO 52 I-I1#LX 

52 X(I) - X( I) +TX 

CALL LINE(X(ID,CP( n)#ND»0»l»C./i.,C.D. ) 
GO TC 21 

61 PEi^IND N1 

M ■ M +1 

CALL PLOT( 12.»0.,--j) 

IF (^.GT.2) GO TC 71 
SX - -SClLX*XC,(KTEi) 

ENC0Ce(2<f»62#R) 

62 FQRMAT( 24HUPPEP SURFACE PRESSURE ) 

CALL SYR<0QL(O., .5,.lA,R,C.,2<t) 
ENCQDE(2<*>6A>R) 

6^ FCRMAT( 2AHL0WER SURFACE PRESSURE ) 

CALL SYM8CL(3.5# .5>.l‘»»R^a.,2^) 

GO TC 11 

71 10 -1 

RETURN 

ICl 10 ■ C 

CALL PL0T(12.,C.,-5) 

RETURN 

FNO 
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I 

ft 

This report was prepared as an account of 
Government sponsored work. Neither the 
United States, nor the Administration, 
nor any person acting on behalf of the 
Administration: 

A. Makes any warranty or representation, 
express or implied, with respect to the 
accuracy, completeness, or usefulness of 
the Information contained in this report, 
or that the use of any information, 
apparatus, method, or process disclosed 
in this report may not infringe privately 
owned rights; or 

B. Assumes any liabilities with respect to 
the use of, or for damages resulting from 
the use of any information, apparatus, 
method, or process disclosed in this 
report . 

As used in the above, "person acting on behalf 
of the Administration" includes any employee 
or contractor of the Administration, or 
employee of such contractor, to the extent 
that such employee or contractor of the 
Administration, or employee of such contractor 
prepares, disseminates, or provides access to, 
any information pursuant to his emplo5nnent or 
contract with the Administration, or his 
employment with such contractor. 








